curcomp <- "~/Dropbox/Boards/boards_rep/"

# load all packages
setwd(curcomp)
suppressMessages(library(pbapply))
suppressMessages(library(parallel))
suppressMessages(library(Matrix))
suppressMessages(library(car))
suppressMessages(library(plm))
suppressMessages(library(xtable))
suppressMessages(library(MASS))
suppressMessages(library(tidyverse))

# Define custom helper functions 
myround <- function(vec, dig){return(formatC(round(as.numeric(as.character(vec)),dig),dig,format="f"))}
ptable <- function(contents, file){
  print(xtable(contents), type='latex', sanitize.text.function=identity, include.rownames = FALSE, include.colnames = FALSE, 
        file = file, only.contents = TRUE, hline.after = NULL, add.to.row = addtorow)
}


######################
# This script conducts the analysis
###### Contents ######
# 1. Load analysis data
# 2. Main models and figure
# 3. Models with firm-level political risk
# 4. Models with lags
# 5. Models with indirect interlocks
# 6. Models with both indirect and direct interlocks
# 7. Models with binary indirect interlocks
# 8. Models with initial adoption only
# 9. Models interacting interlock-weighted variables with sector emissions intensity
# 10. Models with pro-climate lobbying firm interlocks 
# 11. Summary statistics on interlocks
# 12. Main models not dropping no revt/no emp firms

# 1. Load analysis data -------------------------------------
dat <- read_csv("analysis_data_revt_drop.csv")

dat$ind_year <- paste(substr(dat$naics,1,2), dat$year, sep = "")
dat$ind <- substr(dat$naics,1,3)
dat$cso_exists_sewtd <- dat$cso_exists_sewtd*100000
dat$cdp_report_sewtd <- dat$cdp_report_sewtd*100000
dat$numsupcoal_sewtd <- dat$numsupcoal_sewtd*100000
dat$numoppcoal_sewtd <- dat$numoppcoal_sewtd*100000
dat$lob_sewtd <- dat$lob_sewtd*100000
dat$evc <- dat$evc/.001
dat$flpr_envrisk <- dat$flpr_envrisk/max(dat$flpr_envrisk, na.rm = TRUE)
dat$numsupcoal_d <- as.numeric(dat$numsupcoal > 0)
dat$numoppcoal_d <- as.numeric(dat$numoppcoal > 0)
dat$emp <- asinh(dat$emp)
dat$revt <- asinh(dat$revt)

dat$num_interlocks <- asinh(dat$num_interlocks)
dat$cso_exists_ilwtd <- asinh(dat$cso_exists_ilwtd)
dat$cdp_report_ilwtd <- asinh(dat$cdp_report_ilwtd)
dat$numsupcoal_ilwtd <- asinh(dat$numsupcoal_ilwtd)
dat$lob_ilwtd <- asinh(dat$lob_ilwtd)
dat$numoppcoal_ilwtd <- asinh(dat$numoppcoal_ilwtd)

dat$num_indinterlocks <- asinh(dat$num_indinterlocks)
dat$cso_exists_indilwtd <- asinh(dat$cso_exists_indilwtd)
dat$cdp_report_indilwtd <- asinh(dat$cdp_report_indilwtd)
dat$numsupcoal_indilwtd <- asinh(dat$numsupcoal_indilwtd)
dat$lob_indilwtd <- asinh(dat$lob_indilwtd)
dat$numoppcoal_indilwtd <- asinh(dat$numoppcoal_indilwtd)

dat$cso_exists_indilwtd_alt <- asinh(dat$cso_exists_indilwtd_alt)
dat$cdp_report_indilwtd_alt <- asinh(dat$cdp_report_indilwtd_alt)
dat$numsupcoal_indilwtd_alt <- asinh(dat$numsupcoal_indilwtd_alt)
dat$lob_indilwtd_alt <- asinh(dat$lob_indilwtd_alt)
dat$numoppcoal_indilwtd_alt <- asinh(dat$numoppcoal_indilwtd_alt)

# Create pdata.frame and demean by firm 
pdat <- pdata.frame(dat,index = c("gvkey","year"))
varlist <- colnames(pdat)[-which(colnames(pdat) %in% c("gvkey","conml","year","cik","cusip","naics","emissions_direct_intensity","ind_year","ind"))]
for(v in varlist){
  pdat[,v] <- Within(pdat[,v],effect = "individual",na.rm = T)
}

# Subset curtailed data frame
cur <- dat[dat$year %in% 2004:2021,]
# Create pdata.frame and demean by firm 
pcur <- pdata.frame(cur,index = c("gvkey","year"))
varlist <- colnames(pdat)[-which(colnames(pdat) %in% c("gvkey","conml","year","cik","cusip","naics","emissions_direct_intensity","ind_year","ind"))]
for(v in varlist){
  pcur[,v] <- Within(pcur[,v],effect = "individual",na.rm = T)
}

# 2. Main models -------------------------------------
mod1 <- lm(cso_exists ~ cso_exists_ilwtd, data = cur)
mod1aa <- lm(cso_exists ~ cso_exists_ilwtd + cso_exists_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = cur)
mod1a <- lm(cso_exists ~ cso_exists_ilwtd + cso_exists_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = cur)
mod1b <- lm(cso_exists ~ cso_exists_ilwtd + cso_exists_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = pcur)

mod2 <- lm(cdp_report ~ cdp_report_ilwtd, data = dat)
mod2aa <- lm(cdp_report ~ cdp_report_ilwtd + cdp_report_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat)
mod2a <- lm(cdp_report ~ cdp_report_ilwtd + cdp_report_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = dat)
mod2b <- lm(cdp_report ~ cdp_report_ilwtd + cdp_report_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = pdat)

mod3 <- lm(I(numsupcoal>0) ~ numsupcoal_ilwtd, data = dat)
mod3aa <- lm(I(numsupcoal>0) ~ numsupcoal_ilwtd + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat)
mod3a <- lm(I(numsupcoal>0) ~ numsupcoal_ilwtd + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = dat)
mod3b <- lm(numsupcoal_d ~ numsupcoal_ilwtd + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = pdat)

mod4 <- lm(lob ~ lob_ilwtd, data = dat)
mod4aa <- lm(lob ~ lob_ilwtd + lob_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat)
mod4a <- lm(lob ~ lob_ilwtd + lob_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = dat)
mod4b <- lm(lob ~ lob_ilwtd + lob_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = pdat)

mod5 <- lm(I(numoppcoal>0) ~ numoppcoal_ilwtd, data = dat)
mod5aa <- lm(I(numoppcoal>0) ~ numoppcoal_ilwtd + numoppcoal_sewtd + num_interlocks + evc + numsupcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat)
mod5a <- lm(I(numoppcoal>0) ~ numoppcoal_ilwtd + numoppcoal_sewtd + num_interlocks + evc + numsupcoal + emp + revt + as.factor(ind_year), data = dat)
mod5b <- lm(numoppcoal_d ~ numoppcoal_ilwtd + numoppcoal_sewtd + num_interlocks + evc + numsupcoal + emp + revt + as.factor(ind_year), data = pdat)

# GLMs for Table 2
mod1_glm <- glm(cso_exists ~ cso_exists_ilwtd, data = cur, family = "binomial")
mod1aa_glm <- glm(cso_exists ~ cso_exists_ilwtd + cso_exists_sewtd + num_interlocks + numoppcoal + emp + revt, data = cur, family = "binomial")
mod2_glm <- glm(cdp_report ~ cdp_report_ilwtd, data = dat, family = "binomial")
mod2aa_glm <- glm(cdp_report ~ cdp_report_ilwtd + cdp_report_sewtd + num_interlocks + evc + numoppcoal + emp + revt, data = dat, family = "binomial")
mod3_glm <- glm(I(numsupcoal>0) ~ numsupcoal_ilwtd, data = dat, maxit = 100, family = "binomial")
mod3aa_glm <- glm(I(numsupcoal>0) ~ numsupcoal_ilwtd + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + emp + revt, data = dat, maxit = 100, family = "binomial")
mod4_glm <- glm(lob ~ lob_ilwtd, data = dat, family = "binomial")
mod4aa_glm <- glm(lob ~ lob_ilwtd + lob_sewtd + num_interlocks + evc + numoppcoal + emp + revt, data = dat, family = "binomial")
mod5_glm <- glm(I(numoppcoal>0) ~ numoppcoal_ilwtd, data = dat, family = "binomial")
mod5aa_glm <- glm(I(numoppcoal>0) ~ numoppcoal_ilwtd + numoppcoal_sewtd + num_interlocks + evc + numsupcoal + emp + revt, data = dat, family = "binomial")

# Write Table 1 
mods <- list(summary(mod1b)$coef[2:8,],summary(mod2b)$coef[2:8,],summary(mod3b)$coef[2:8,],summary(mod4b)$coef[2:8,])
rownames(mods[[1]])[2] <- rownames(mods[[2]])[2] <- rownames(mods[[3]])[2] <- rownames(mods[[4]])[2] <- "Structural"
vars <- c(rownames(mods[[1]])[1],rownames(mods[[2]])[1],rownames(mods[[3]])[1],rownames(mods[[4]])[1],rownames(mods[[1]])[2:7])
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17,19)] <- 
  c("CSOs, Interlock wtd.","CDP reporting, Interlock wtd.","Climate coalitions, Interlock wtd.","Climate lobbying, Interlock wtd.",
    "Structural eqv. wtd. DV","Num. interlocks", "Eig. Centrality", "Num. opp. coalitions","Employees","Revenue") 
N <- c(sum(summary(mod1b)$df[1:2]),sum(summary(mod2b)$df[1:2]),sum(summary(mod3b)$df[1:2]),sum(summary(mod4b)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(20,21); addtorow$command <- c(paste0(' '),paste0(''))
ptable(cbind(rownames(tab), tab), "output/main_saturated.tex")

# Table 2, A1
coefs_CSO <- c(summary(mod1)$coef[2,1],summary(mod1aa)$coef[2,1],summary(mod1a)$coef[2,1],summary(mod1b)$coef[2,1],summary(mod1_glm)$coef[2,1],summary(mod1aa_glm)$coef[2,1])
se_CSO <- c(summary(mod1)$coef[2,2],summary(mod1aa)$coef[2,2],summary(mod1a)$coef[2,2],summary(mod1b)$coef[2,2],summary(mod1_glm)$coef[2,2],summary(mod1aa_glm)$coef[2,2])
coefs_CSO[1:4] <- 100*coefs_CSO[1:4]; se_CSO[1:4] <- 100*se_CSO[1:4]
se_CSO <- paste("(", myround(se_CSO,2), ")", sep = "")
p_CSO <- c(summary(mod1)$coef[2,4],summary(mod1aa)$coef[2,4],summary(mod1a)$coef[2,4],summary(mod1b)$coef[2,4],summary(mod1_glm)$coef[2,4],summary(mod1aa_glm)$coef[2,4])
stars <- rep("", length(coefs_CSO)); stars[p_CSO < .1] <- "^{+}"; stars[p_CSO < .05] <- "^{*}"; stars[p_CSO < .01] <- "^{**}"; stars[p_CSO < .001] <- "^{***}";  
coefs_CSO <- paste(myround(coefs_CSO,2), stars, sep = "")

coefs_CDP <- c(summary(mod2)$coef[2,1],summary(mod2aa)$coef[2,1],summary(mod2a)$coef[2,1],summary(mod2b)$coef[2,1],summary(mod2_glm)$coef[2,1],summary(mod2aa_glm)$coef[2,1])
se_CDP <- c(summary(mod2)$coef[2,2],summary(mod2aa)$coef[2,2],summary(mod2a)$coef[2,2],summary(mod2b)$coef[2,2],summary(mod2_glm)$coef[2,2],summary(mod2aa_glm)$coef[2,2])
coefs_CDP[1:4] <- 100*coefs_CDP[1:4]; se_CDP[1:4] <- 100*se_CDP[1:4]
se_CDP <- paste("(", myround(se_CDP,2), ")", sep = "")
p_CDP <- c(summary(mod2)$coef[2,4],summary(mod2aa)$coef[2,4],summary(mod2a)$coef[2,4],summary(mod2b)$coef[2,4],summary(mod2_glm)$coef[2,4],summary(mod2aa_glm)$coef[2,4])
stars <- rep("", length(coefs_CDP)); stars[p_CDP < .1] <- "^{+}"; stars[p_CDP < .05] <- "^{*}"; stars[p_CDP < .01] <- "^{**}"; stars[p_CDP < .001] <- "^{***}";  
coefs_CDP <- paste(myround(coefs_CDP,2), stars, sep = "")

coefs_COL <- c(summary(mod3)$coef[2,1],summary(mod3aa)$coef[2,1],summary(mod3a)$coef[2,1],summary(mod3b)$coef[2,1],summary(mod3_glm)$coef[2,1],summary(mod3aa_glm)$coef[2,1])
se_COL <- c(summary(mod3)$coef[2,2],summary(mod3aa)$coef[2,2],summary(mod3a)$coef[2,2],summary(mod3b)$coef[2,2],summary(mod3_glm)$coef[2,2],summary(mod3aa_glm)$coef[2,2])
coefs_COL[1:4] <- 100*coefs_COL[1:4]; se_COL[1:4] <- 100*se_COL[1:4]
se_COL <- paste("(", myround(se_COL,2), ")", sep = "")
p_COL <- c(summary(mod3)$coef[2,4],summary(mod3aa)$coef[2,4],summary(mod3a)$coef[2,4],summary(mod3b)$coef[2,4],summary(mod3_glm)$coef[2,4],summary(mod3aa_glm)$coef[2,4])
stars <- rep("", length(coefs_COL)); stars[p_COL < .1] <- "^{+}"; stars[p_COL < .05] <- "^{*}"; stars[p_COL < .01] <- "^{**}"; stars[p_COL < .001] <- "^{***}";  
coefs_COL <- paste(myround(coefs_COL,2), stars, sep = "")

coefs_LOB <- c(summary(mod4)$coef[2,1],summary(mod4aa)$coef[2,1],summary(mod4a)$coef[2,1],summary(mod4b)$coef[2,1],summary(mod4_glm)$coef[2,1],summary(mod4aa_glm)$coef[2,1])
se_LOB <- c(summary(mod4)$coef[2,2],summary(mod4aa)$coef[2,2],summary(mod4a)$coef[2,2],summary(mod4b)$coef[2,2],summary(mod4_glm)$coef[2,2],summary(mod4aa_glm)$coef[2,2])
coefs_LOB[1:4] <- 100*coefs_LOB[1:4]; se_LOB[1:4] <- 100*se_LOB[1:4]
se_LOB <- paste("(", myround(se_LOB,2), ")", sep = "")
p_LOB <- c(summary(mod4)$coef[2,4],summary(mod4aa)$coef[2,4],summary(mod4a)$coef[2,4],summary(mod4b)$coef[2,4],summary(mod4_glm)$coef[2,4],summary(mod4aa_glm)$coef[2,4])
stars <- rep("", length(coefs_LOB)); stars[p_LOB < .1] <- "^{+}"; stars[p_LOB < .05] <- "^{*}"; stars[p_LOB < .01] <- "^{**}"; stars[p_LOB < .001] <- "^{***}";  
coefs_LOB <- paste(myround(coefs_LOB,2), stars, sep = "")

coefs_OPP <- c(summary(mod5)$coef[2,1],summary(mod5aa)$coef[2,1],summary(mod5a)$coef[2,1],summary(mod5b)$coef[2,1],summary(mod5_glm)$coef[2,1],summary(mod5aa_glm)$coef[2,1])
se_OPP <- c(summary(mod5)$coef[2,2],summary(mod5aa)$coef[2,2],summary(mod5a)$coef[2,2],summary(mod5b)$coef[2,2],summary(mod5_glm)$coef[2,2],summary(mod5aa_glm)$coef[2,2])
coefs_OPP[1:4] <- 100*coefs_OPP[1:4]; se_OPP[1:4] <- 100*se_OPP[1:4]
se_OPP <- paste("(", myround(se_OPP,2), ")", sep = "")
p_OPP <- c(summary(mod5)$coef[2,4],summary(mod5aa)$coef[2,4],summary(mod5a)$coef[2,4],summary(mod5b)$coef[2,4],summary(mod5_glm)$coef[2,4],summary(mod5aa_glm)$coef[2,4])
stars <- rep("", length(coefs_OPP)); stars[p_OPP < .1] <- "^{+}"; stars[p_OPP < .05] <- "^{*}"; stars[p_OPP < .01] <- "^{**}"; stars[p_OPP < .001] <- "^{***}";  
coefs_OPP <- paste(myround(coefs_OPP,2), stars, sep = "")

addtorow <- list(); addtorow$pos <- list(1); addtorow$command <- c("")
ptable(rbind(c("CSO, Interlock wtd.", coefs_CSO),c("", se_CSO)), "output/tab2_CSO.tex")
addtorow <- list(); addtorow$pos <- list(1); addtorow$command <- c("")
ptable(rbind(c("CDP reporting, Interlock wtd.", coefs_CDP),c("", se_CDP)), "output/tab2_CDP.tex")
addtorow <- list(); addtorow$pos <- list(1); addtorow$command <- c("")
ptable(rbind(c("Climate coalitions, Interlock wtd.", coefs_COL),c("", se_COL)), "output/tab2_COL.tex")
addtorow <- list(); addtorow$pos <- list(1); addtorow$command <- c("")
ptable(rbind(c("Climate lobbying, Interlock wtd.", coefs_LOB),c("", se_LOB)), "output/tab2_LOB.tex")
addtorow <- list(); addtorow$pos <- list(1); addtorow$command <- c("")
ptable(rbind(c("Opposing coalitions, Interlock wtd.", coefs_OPP),c("", se_OPP)), "output/tab2_OPP.tex")

# clean up
rems <- ls()[ls() %in% c("dat","cur","pdat","pcur","myround","ptable","dat_lag","cur_lag") == FALSE]
rm(list = rems)

## Figure 1
set.seed(12345)
quantile(cur$cso_exists_ilwtd[cur$cso_exists_ilwtd > 0], c(.5,.8,.9), na.rm = TRUE)

pcur_illus <- pcur
pcur_illus$cso_exists_sewtd <- pcur_illus$cso_exists_sewtd/100000
mod1b_illus <- lm(cso_exists ~ cso_exists_ilwtd + cso_exists_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), 
  data = pcur_illus)

# cso_exists figure.
ests <- mod1b_illus$coef; ests_cv <- vcov(mod1b_illus); 
betas <- mvrnorm(1000, ests, ests_cv)
cur1 <- cur0 <- as.matrix(model.matrix(mod1b_illus))
cur0[,"cso_exists_ilwtd"] <- 0; cur1[,"cso_exists_ilwtd"] <- asinh(1)
agg <- aggregate(cso_exists ~ as.character(gvkey), data = cur, FUN = mean, na.rm = TRUE); colnames(agg) <- c("gvkey","cso_exists_firmavg")
cur0 <- data.frame(cur0)
cur0$gvkey <- c(unlist(strsplit(rownames(cur0), "-")))[2*(1:nrow(cur0)) - 1]
cur0 <- left_join(data.frame(cur0), agg)
cso_exists_firmavg <- cur0$cso_exists_firmavg
cur0 <- as.matrix(cur0[,-c(ncol(cur0)-1,ncol(cur0))])
preds0_cso <- 100*(apply(cur0%*%t(betas), 2, mean) + mean(cso_exists_firmavg)); 
preds1_cso <- 100*(apply(cur1%*%t(betas), 2, mean) + mean(cso_exists_firmavg)); 
mean(preds0_cso)
mean(preds1_cso)
mean(preds1_sup) - mean(preds0_sup)
quantile(cur$cdp_report_ilwtd[cur$cdp_report_ilwtd > 0], c(.5,.8,.9), na.rm = TRUE)

pdat_illus <- pdat
# pdat_illus$cdp_report_sewtd <- pdat_illus$cdp_report_sewtd*1e16
mod2b_illus <- lm(cdp_report ~ cdp_report_ilwtd + cdp_report_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), 
                  data = pdat_illus)

# cdp_report figure.
ests <- mod2b_illus$coef; ests_cv <- vcov(mod2b_illus); 
betas <- mvrnorm(1000, ests, ests_cv)
cur1 <- cur0 <- model.matrix(mod2b_illus)
cur0[,"cdp_report_ilwtd"] <- 0; cur1[,"cdp_report_ilwtd"] <- 1
agg <- aggregate(cdp_report ~ as.character(gvkey), data = dat, FUN = mean, na.rm = TRUE); colnames(agg) <- c("gvkey","cdp_report_firmavg")
cur0 <- data.frame(cur0)
cur0$gvkey <- c(unlist(strsplit(rownames(cur0), "-")))[2*(1:nrow(cur0)) - 1]
cur0 <- left_join(data.frame(cur0), agg)
cdp_report_firmavg <- cur0$cdp_report_firmavg
cur0 <- as.matrix(cur0[,-c(ncol(cur0)-1,ncol(cur0))])
preds0_cdp <- 100*(apply(cur0%*%t(betas), 2, mean) + mean(cdp_report_firmavg)); 
preds1_cdp <- 100*(apply(cur1%*%t(betas), 2, mean) + mean(cdp_report_firmavg)); 
mean(preds0_cdp)
mean(preds1_cdp)
mean(preds1_sup) - mean(preds0_sup)

quantile(cur$numsupcoal_ilwtd[cur$numsupcoal_ilwtd > 0], c(.5,.8,.9), na.rm = TRUE)

# pdat_illus$numsupcoal_sewtd <- pdat_illus$numsupcoal_sewtd*1e15
mod3b_illus <- lm(numsupcoal_d ~ numsupcoal_ilwtd + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), 
                  data = pdat_illus)

# numsupcoal figure.
ests <- mod3b_illus$coef; ests_cv <- vcov(mod3b_illus); 
betas <- mvrnorm(1000, ests, ests_cv)
cur1 <- cur0 <- model.matrix(mod3b_illus)
cur0[,"numsupcoal_ilwtd"] <- 0; cur1[,"numsupcoal_ilwtd"] <- 1
agg <- aggregate(I(numsupcoal>0) ~ as.character(gvkey), data = dat, FUN = mean, na.rm = TRUE); colnames(agg) <- c("gvkey","numsupcoal_firmavg")
cur0 <- data.frame(cur0)
cur0$gvkey <- c(unlist(strsplit(rownames(cur0), "-")))[2*(1:nrow(cur0)) - 1]
cur0 <- left_join(data.frame(cur0), agg)
numsupcoal_firmavg <- cur0$numsupcoal_firmavg
cur0 <- as.matrix(cur0[,-c(ncol(cur0)-1,ncol(cur0))])
preds0_sup <- 100*(apply(cur0%*%t(betas), 2, mean) + mean(numsupcoal_firmavg)); 
preds1_sup <- 100*(apply(cur1%*%t(betas), 2, mean) + mean(numsupcoal_firmavg)); 
mean(preds0_sup)
mean(preds1_sup)
mean(preds1_sup) - mean(preds0_sup)

quantile(cur$numsupcoal_ilwtd[cur$numsupcoal_ilwtd > 0], c(.5,.8,.9), na.rm = TRUE)

# pdat_illus$lob_sewtd <- pdat_illus$lob_sewtd*1e16
mod4b_illus <- lm(lob ~ lob_ilwtd + lob_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), 
                  data = pdat_illus)

# lob figure.
ests <- mod4b_illus$coef; ests_cv <- vcov(mod4b_illus); 
betas <- mvrnorm(1000, ests, ests_cv)
cur1 <- cur0 <- model.matrix(mod4b_illus)
cur0[,"lob_ilwtd"] <- 0; cur1[,"lob_ilwtd"] <- 1
agg <- aggregate(lob ~ as.character(gvkey), data = dat, FUN = mean, na.rm = TRUE); colnames(agg) <- c("gvkey","lob_firmavg")
cur0 <- data.frame(cur0)
cur0$gvkey <- c(unlist(strsplit(rownames(cur0), "-")))[2*(1:nrow(cur0)) - 1]
cur0 <- left_join(data.frame(cur0), agg)
lob_firmavg <- cur0$lob_firmavg
cur0 <- as.matrix(cur0[,-c(ncol(cur0)-1,ncol(cur0))])
preds0_lob <- 100*(apply(cur0%*%t(betas), 2, mean) + mean(lob_firmavg)); 
preds1_lob <- 100*(apply(cur1%*%t(betas), 2, mean) + mean(lob_firmavg)); 
mean(preds0_lob)
mean(preds1_lob)
mean(preds1_lob) - mean(preds0_lob)

quantile(cur$lob_ilwtd[cur$lob_ilwtd > 0], c(.5,.8,.9), na.rm = TRUE)

pdf("output/table1_illustrated.pdf", height = 3.8, width = 6)
set.seed(12345)
par(mfrow = c(1,1), mai = c(.1, 0.3, .11, .3), mar  = c(1, 2, 1, 0))
plot(-5,-5, xlim = c(-.02,4.3), ylim = c(-.3,7.8), xlab = "", ylab = "", axes = F)
axis(2, cex.axis = .6, mgp=c(3, .45, 0), tck=-0.02, labels = as.character(0:7), at = 0:7)
title(main = "Predicted outcomes increasing interlock wtd. measures from 0 to 1", cex.main = .65, line = .1)
title(ylab = "Percentage chance of undertaking activity", cex.lab = .65, line = 1.2)
rect(xleft = .2, xright = .4, ybottom = 0, ytop = mean(preds0_cso))
rect(xleft = .6, xright = .8, ybottom = 0, ytop = mean(preds1_cso))
ys0 <- quantile(preds0_cso + rnorm(1000, 0, summary(mod1b_illus)$sigma), c(.025,.975))
segments(x0 = .3, x1 = .3, y0 = ys0[1], y1 = ys0[2])
ys1 <- quantile(preds1_cso + rnorm(1000, 0, summary(mod1b_illus)$sigma), c(.025,.975))
segments(x0 = .7, x1 = .7, y0 = ys1[1], y1 = ys1[2])
text(x = c(.3,.7), y = -.2, labels = c("0","1"), cex = .5)
text(x = c(.5), y = -.45, labels = c("Interlock Wtd. CSOs"), cex = .4)
text(x = .5, y = 7.3, labels = c("CSO"), cex = .6)

quantile((preds1_cso - preds0_cso)/preds0_cso)

rect(xleft = 1.2, xright = 1.4, ybottom = 0, ytop = mean(preds0_cdp))
rect(xleft = 1.6, xright = 1.8, ybottom = 0, ytop = mean(preds1_cdp))
ys0 <- quantile(preds0_cdp + rnorm(1000, 0, summary(mod2b_illus)$sigma), c(.025,.975))
segments(x0 = 1.3, x1 = 1.3, y0 = ys0[1], y1 = ys0[2])
ys1 <- quantile(preds1_cdp + rnorm(1000, 0, summary(mod2b_illus)$sigma), c(.025,.975))
segments(x0 = 1.7, x1 = 1.7, y0 = ys1[1], y1 = ys1[2])
text(x = c(1.3,1.7), y = -.2, labels = c("0","1"), cex = .5)
text(x = c(1.5), y = -.45, labels = c("Interlock Wtd. CDP Reporting"), cex = .4)
text(x = 1.5, y = 7.3, labels = c("CDP reporting"), cex = .6)

rect(xleft = 2.2, xright = 2.4, ybottom = 0, ytop = mean(preds0_sup))
rect(xleft = 2.6, xright = 2.8, ybottom = 0, ytop = mean(preds1_sup))
ys0 <- quantile(preds0_sup + rnorm(1000, 0, summary(mod3b_illus)$sigma), c(.025,.975))
segments(x0 = 2.3, x1 = 2.3, y0 = ys0[1], y1 = ys0[2])
ys1 <- quantile(preds1_sup + rnorm(1000, 0, summary(mod3b_illus)$sigma), c(.025,.975))
segments(x0 = 2.7, x1 = 2.7, y0 = ys1[1], y1 = ys1[2])
text(x = c(2.3,2.7), y = -.2, labels = c("0","1"), cex = .5)
text(x = c(2.5), y = -.45, labels = c("Interlock Wtd. Coalitions"), cex = .4)
text(x = 2.5, y = 7.3, labels = c("Climate coalitions"), cex = .6)

rect(xleft = 3.2, xright = 3.4, ybottom = 0, ytop = mean(preds0_lob))
rect(xleft = 3.6, xright = 3.8, ybottom = 0, ytop = mean(preds1_lob))
ys0 <- quantile(preds0_lob + rnorm(1000, 0, summary(mod4b_illus)$sigma), c(.025,.975))
segments(x0 = 3.3, x1 = 3.3, y0 = ys0[1], y1 = ys0[2])
ys1 <- quantile(preds1_lob + rnorm(1000, 0, summary(mod4b_illus)$sigma), c(.025,.975))
segments(x0 = 3.7, x1 = 3.7, y0 = ys1[1], y1 = ys1[2])
text(x = c(3.3,3.7), y = -.2, labels = c("0","1"), cex = .5)
text(x = c(3.5), y = -.45, labels = c("Interlock Wtd. Lobbying"), cex = .4)
text(x = 3.5, y = 7.3, labels = c("Lobbying"), cex = .6)
dev.off()

# clean up
rems <- ls()[ls() %in% c("dat","cur","pdat","pcur","myround","ptable","dat_lag","cur_lag") == FALSE]
rm(list = rems)

# 3. Models with firm-level political risk -------------------------------------
mod1aa_flpr <- lm(cso_exists ~ cso_exists_ilwtd + cso_exists_sewtd + num_interlocks + evc + flpr_envrisk + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = cur)
mod1a_flpr <- lm(cso_exists ~ cso_exists_ilwtd + cso_exists_sewtd + num_interlocks + evc + flpr_envrisk + numoppcoal + emp + revt + as.factor(ind_year), data = cur)
mod1b_flpr <- lm(cso_exists ~ cso_exists_ilwtd + cso_exists_sewtd + num_interlocks + evc + flpr_envrisk + numoppcoal + emp + revt + as.factor(ind_year), data = pcur)

mod2aa_flpr <- lm(cdp_report ~ cdp_report_ilwtd + cdp_report_sewtd + num_interlocks + evc + flpr_envrisk + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat)
mod2a_flpr <- lm(cdp_report ~ cdp_report_ilwtd + cdp_report_sewtd + num_interlocks + evc + flpr_envrisk + numoppcoal + emp + revt + as.factor(ind_year), data = dat)
mod2b_flpr <- lm(cdp_report ~ cdp_report_ilwtd + cdp_report_sewtd + num_interlocks + evc + flpr_envrisk + numoppcoal + emp + revt + as.factor(ind_year), data = pdat)

mod3aa_flpr <- lm(I(numsupcoal>0) ~ numsupcoal_ilwtd + numsupcoal_sewtd + num_interlocks + evc + flpr_envrisk + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat)
mod3a_flpr <- lm(I(numsupcoal>0) ~ numsupcoal_ilwtd + numsupcoal_sewtd + num_interlocks + evc + flpr_envrisk + numoppcoal + emp + revt + as.factor(ind_year), data = dat)
mod3b_flpr <- lm(numsupcoal_d ~ numsupcoal_ilwtd + numsupcoal_sewtd + num_interlocks + evc + flpr_envrisk + numoppcoal + emp + revt + as.factor(ind_year), data = pdat)

mod4aa_flpr <- lm(lob ~ lob_ilwtd + lob_sewtd + num_interlocks + evc + flpr_envrisk + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat)
mod4a_flpr <- lm(lob ~ lob_ilwtd + lob_sewtd + num_interlocks + evc + flpr_envrisk + numoppcoal + emp + revt + as.factor(ind_year), data = dat)
mod4b_flpr <- lm(lob ~ lob_ilwtd + lob_sewtd + num_interlocks + evc + flpr_envrisk + numoppcoal + emp + revt + as.factor(ind_year), data = pdat)

mods <- list(summary(mod1aa_flpr)$coef[2:9,],summary(mod2aa_flpr)$coef[2:9,],summary(mod3aa_flpr)$coef[2:9,],summary(mod4aa_flpr)$coef[2:9,])
rownames(mods[[1]])[2] <- rownames(mods[[2]])[2] <- rownames(mods[[3]])[2] <- rownames(mods[[4]])[2] <- "Structural"
vars <- c(rownames(mods[[1]])[1],rownames(mods[[2]])[1],rownames(mods[[3]])[1],rownames(mods[[4]])[1],rownames(mods[[1]])[2:8])
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17,19,21)] <- 
  c("CSOs, Interlock wtd.","CDP reporting, Interlock wtd.","Climate coalitions, Interlock wtd.","Climate lobbying, Interlock wtd.",
    "Structural eqv. wtd. DV","Num. interlocks", "Eig. Centrality","Env. risk", "Num. opp. coalitions","Employees","Revenue") 
N <- c(sum(summary(mod1aa_flpr)$df[1:2]),sum(summary(mod2aa_flpr)$df[1:2]),sum(summary(mod3aa_flpr)$df[1:2]),sum(summary(mod4aa_flpr)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(22,23); addtorow$command <- c(paste0(' '),paste0(''))
ptable(cbind(rownames(tab), tab), "output/flpr_indnotyearfe.tex")

mods <- list(summary(mod1a_flpr)$coef[2:9,],summary(mod2a_flpr)$coef[2:9,],summary(mod3a_flpr)$coef[2:9,],summary(mod4a_flpr)$coef[2:9,])
rownames(mods[[1]])[2] <- rownames(mods[[2]])[2] <- rownames(mods[[3]])[2] <- rownames(mods[[4]])[2] <- "Structural"
vars <- c(rownames(mods[[1]])[1],rownames(mods[[2]])[1],rownames(mods[[3]])[1],rownames(mods[[4]])[1],rownames(mods[[1]])[2:8])
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17,19,21)] <- 
  c("CSOs, Interlock wtd.","CDP reporting, Interlock wtd.","Climate coalitions, Interlock wtd.","Climate lobbying, Interlock wtd.",
    "Structural eqv. wtd. DV","Num. interlocks", "Eig. Centrality","Env. risk", "Num. opp. coalitions","Employees","Revenue") 
N <- c(sum(summary(mod1a_flpr)$df[1:2]),sum(summary(mod2a_flpr)$df[1:2]),sum(summary(mod3a_flpr)$df[1:2]),sum(summary(mod4a_flpr)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(22,23); addtorow$command <- c(paste0(' '),paste0(''))
ptable(cbind(rownames(tab), tab), "output/flpr_indfe.tex")

mods <- list(summary(mod1b_flpr)$coef[2:9,],summary(mod2b_flpr)$coef[2:9,],summary(mod3b_flpr)$coef[2:9,],summary(mod4b_flpr)$coef[2:9,])
rownames(mods[[1]])[2] <- rownames(mods[[2]])[2] <- rownames(mods[[3]])[2] <- rownames(mods[[4]])[2] <- "Structural"
vars <- c(rownames(mods[[1]])[1],rownames(mods[[2]])[1],rownames(mods[[3]])[1],rownames(mods[[4]])[1],rownames(mods[[1]])[2:8])
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17,19,21)] <- 
  c("CSOs, Interlock wtd.","CDP reporting, Interlock wtd.","Climate coalitions, Interlock wtd.","Climate lobbying, Interlock wtd.",
    "Structural eqv. wtd. DV","Num. interlocks", "Eig. Centrality","Env. risk", "Num. opp. coalitions","Employees","Revenue") 
N <- c(sum(summary(mod1b_flpr)$df[1:2]),sum(summary(mod2b_flpr)$df[1:2]),sum(summary(mod3b_flpr)$df[1:2]),sum(summary(mod4b_flpr)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(22,23); addtorow$command <- c(paste0(' '),paste0(''))
ptable(cbind(rownames(tab), tab), "output/flpr_saturated.tex")

# clean up
rems <- ls()[ls() %in% c("dat","cur","pdat","pcur","myround","ptable","dat_lag","cur_lag") == FALSE]
rm(list = rems)

# 4. Models with lags -------------------------------------
dat_lag <- dat %>%
  group_by(gvkey) %>%
  mutate(year = year,
         cso_exists_lag = lag(cso_exists,1),cso_exists_ilwtd_lag = lag(cso_exists_ilwtd,1),cso_exists_sewtd_lag = lag(cso_exists_sewtd,1),
         cdp_report_lag = lag(cdp_report,1),cdp_report_ilwtd_lag = lag(cdp_report_ilwtd,1),cdp_report_sewtd_lag = lag(cdp_report_sewtd,1),
         numsupcoal_lag = lag(numsupcoal,1),numsupcoal_d_lag = lag(numsupcoal,1),numsupcoal_ilwtd_lag = lag(numsupcoal_ilwtd,1),numsupcoal_sewtd_lag = lag(numsupcoal_sewtd,1),
         lob_lag = lag(lob,1),lob_ilwtd_lag = lag(lob_ilwtd,1),lob_sewtd_lag = lag(lob_sewtd,1))

# Create pdata.frame and demean by firm 
pdat_lag <- pdata.frame(dat_lag,index = c("gvkey","year"))
varlist <- colnames(pdat)[-which(colnames(pdat) %in% c("gvkey","conml","year","cik","cusip","naics","emissions_direct_intensity", "ind_year","ind"))]
pdat_lag <- pdat_lag %>% ungroup()
for(v in varlist){
  pdat_lag[,v] <- Within(pdat_lag[,v],effect = "individual",na.rm = T)
}

# Subset curtailed data frame
cur_lag <- dat_lag[dat_lag$year %in% 2004:2021,]
# Create pdata.frame and demean by firm 
pcur_lag <- pdata.frame(cur_lag,index = c("gvkey","year"))
varlist <- colnames(pdat)[-which(colnames(pdat) %in% c("gvkey","conml","year","cik","cusip","naics","emissions_direct_intensity", "ind_year","ind"))]
for(v in varlist){
  pcur_lag[,v] <- Within(pcur_lag[,v],effect = "individual",na.rm = T)
}

mod1_lag <- lm(cso_exists ~ cso_exists_ilwtd + cso_exists_lag, data = cur_lag)
mod1aa_lag <- lm(cso_exists ~ cso_exists_ilwtd + cso_exists_lag + cso_exists_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = cur_lag)
mod1a_lag <- lm(cso_exists ~ cso_exists_ilwtd + cso_exists_lag + cso_exists_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = cur_lag)
mod1b_lag <- lm(cso_exists ~ cso_exists_ilwtd + cso_exists_lag + cso_exists_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = pcur_lag)

mod2_lag <- lm(cdp_report ~ cdp_report_ilwtd + cdp_report_lag, data = dat_lag)
mod2aa_lag <- lm(cdp_report ~ cdp_report_ilwtd + cdp_report_lag + cdp_report_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat_lag)
mod2a_lag <- lm(cdp_report ~ cdp_report_ilwtd + cdp_report_lag + cdp_report_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = dat_lag)
mod2b_lag <- lm(cdp_report ~ cdp_report_ilwtd + cdp_report_lag + cdp_report_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = pdat_lag)

mod3_lag <- lm(I(numsupcoal>0) ~ numsupcoal_ilwtd + numsupcoal_lag + numsupcoal_sewtd, data = dat_lag)
mod3aa_lag <- lm(I(numsupcoal>0) ~ numsupcoal_ilwtd + numsupcoal_lag + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat_lag)
mod3a_lag <- lm(I(numsupcoal>0) ~ numsupcoal_ilwtd + numsupcoal_lag + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = dat_lag)
mod3b_lag <- lm(numsupcoal_d ~ numsupcoal_ilwtd + numsupcoal_d_lag + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = pdat_lag)

mod4_lag <- lm(lob ~ lob_ilwtd + lob_lag + lob_sewtd, data = dat_lag)
mod4aa_lag <- lm(lob ~ lob_ilwtd + lob_lag + lob_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat_lag)
mod4a_lag <- lm(lob ~ lob_ilwtd + lob_lag + lob_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = dat_lag)
mod4b_lag <- lm(lob ~ lob_ilwtd + lob_lag + lob_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = pdat_lag)

mod1_fd <- lm(I(cso_exists-cso_exists_lag) ~ I(cso_exists_ilwtd-cso_exists_ilwtd_lag), data = cur_lag)
mod2_fd <- lm(I(cdp_report- cdp_report_lag) ~ I(cdp_report_ilwtd-cdp_report_ilwtd_lag), data = dat_lag)
mod3_fd <- lm(I(numsupcoal - numsupcoal_lag) ~ I(numsupcoal_ilwtd - numsupcoal_ilwtd_lag), data = dat_lag)
mod4_fd <- lm(I(lob - lob_lag)~ I(lob_ilwtd - lob_ilwtd_lag), data = dat_lag)

mod1a_fd <- lm(I(cso_exists-cso_exists_lag) ~ I(cso_exists_ilwtd-cso_exists_ilwtd_lag) + cso_exists_sewtd + num_interlocks + evc + numoppcoal + emp + revt, data = cur_lag)
mod2a_fd <- lm(I(cdp_report- cdp_report_lag) ~ I(cdp_report_ilwtd-cdp_report_ilwtd_lag) + cdp_report_sewtd + num_interlocks + evc + numoppcoal + emp + revt, data = dat_lag)
mod3a_fd <- lm(I(numsupcoal - numsupcoal_lag) ~ I(numsupcoal_ilwtd - numsupcoal_ilwtd_lag) + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + emp + revt, data = dat_lag)
mod4a_fd <- lm(I(lob - lob_lag)~ I(lob_ilwtd - lob_ilwtd_lag) + lob_sewtd + num_interlocks + evc + numoppcoal + emp + revt, data = dat_lag)

mods <- list(summary(mod1_lag)$coef[1:3,],summary(mod2_lag)$coef[1:3,],summary(mod3_lag)$coef[1:3,],summary(mod4_lag)$coef[1:3,])
rownames(mods[[1]])[3] <- rownames(mods[[2]])[3] <- rownames(mods[[3]])[3] <- rownames(mods[[4]])[3] <- "Lagged DV"
vars <- c(rownames(mods[[1]])[2],rownames(mods[[2]])[2],rownames(mods[[3]])[2],rownames(mods[[4]])[2],"Lagged DV","(Intercept)")
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11)] <- 
  c("CSOs, Interlock wtd.","CDP reporting, Interlock wtd.","Climate coalitions, Interlock wtd.","Climate lobbying, Interlock wtd.","Lagged DV","Intercept") 
N <- c(sum(summary(mod1_lag)$df[1:2]),sum(summary(mod2_lag)$df[1:2]),sum(summary(mod3_lag)$df[1:2]),sum(summary(mod4_lag)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(12,13); addtorow$command <- c(paste0(' '),paste0(' '))
ptable(cbind(rownames(tab), tab), "output/lag_bivariate.tex")

mods <- list(summary(mod1aa_lag)$coef[2:9,],summary(mod2aa_lag)$coef[2:9,],summary(mod3aa_lag)$coef[2:9,],summary(mod4aa_lag)$coef[2:9,])
rownames(mods[[1]])[2] <- rownames(mods[[2]])[2] <- rownames(mods[[3]])[2] <- rownames(mods[[4]])[2] <- "Lagged DV"
rownames(mods[[1]])[3] <- rownames(mods[[2]])[3] <- rownames(mods[[3]])[3] <- rownames(mods[[4]])[3] <- "Structural"
vars <- c(rownames(mods[[1]])[1],rownames(mods[[2]])[1],rownames(mods[[3]])[1],rownames(mods[[4]])[1],rownames(mods[[1]])[2:8])
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17,19,21)] <- 
  c("CSOs, Interlock wtd.","CDP reporting, Interlock wtd.","Climate coalitions, Interlock wtd.","Climate lobbying, Interlock wtd.",
    "Lagged DV", "Structural eqv. wtd. DV","Num. interlocks", "Eig. Centrality", "Num. opp. coalitions","Employees","Revenue") 
N <- c(sum(summary(mod1aa_lag)$df[1:2]),sum(summary(mod2aa_lag)$df[1:2]),sum(summary(mod3aa_lag)$df[1:2]),sum(summary(mod4aa_lag)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(22,23); addtorow$command <- c(paste0(' '),paste0(''))
ptable(cbind(rownames(tab), tab), "output/lag_indnotyearfe.tex")

mods <- list(summary(mod1a_lag)$coef[2:9,],summary(mod2a_lag)$coef[2:9,],summary(mod3a_lag)$coef[2:9,],summary(mod4a_lag)$coef[2:9,])
rownames(mods[[1]])[2] <- rownames(mods[[2]])[2] <- rownames(mods[[3]])[2] <- rownames(mods[[4]])[2] <- "Lagged DV"
rownames(mods[[1]])[3] <- rownames(mods[[2]])[3] <- rownames(mods[[3]])[3] <- rownames(mods[[4]])[3] <- "Structural"
vars <- c(rownames(mods[[1]])[1],rownames(mods[[2]])[1],rownames(mods[[3]])[1],rownames(mods[[4]])[1],rownames(mods[[1]])[2:8])
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17,19,21)] <- 
  c("CSOs, Interlock wtd.","CDP reporting, Interlock wtd.","Climate coalitions, Interlock wtd.","Climate lobbying, Interlock wtd.",
    "Lagged DV","Structural eqv. wtd. DV","Num. interlocks", "Eig. Centrality", "Num. opp. coalitions","Employees","Revenue") 
N <- c(sum(summary(mod1a_lag)$df[1:2]),sum(summary(mod2a_lag)$df[1:2]),sum(summary(mod3a_lag)$df[1:2]),sum(summary(mod4a_lag)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(22,23); addtorow$command <- c(paste0(' '),paste0(''))
ptable(cbind(rownames(tab), tab), "output/lag_indfe.tex")

mods <- list(summary(mod1b_lag)$coef[2:9,],summary(mod2b_lag)$coef[2:9,],summary(mod3b_lag)$coef[2:9,],summary(mod4b_lag)$coef[2:9,])
rownames(mods[[1]])[2] <- rownames(mods[[2]])[2] <- rownames(mods[[3]])[2] <- rownames(mods[[4]])[2] <- "Lagged DV"
rownames(mods[[1]])[3] <- rownames(mods[[2]])[3] <- rownames(mods[[3]])[3] <- rownames(mods[[4]])[3] <- "Structural"
vars <- c(rownames(mods[[1]])[1],rownames(mods[[2]])[1],rownames(mods[[3]])[1],rownames(mods[[4]])[1],rownames(mods[[1]])[2:8])
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17,19,21)] <- 
  c("CSOs, Interlock wtd.","CDP reporting, Interlock wtd.","Climate coalitions, Interlock wtd.","Climate lobbying, Interlock wtd.",
    "Lagged DV","Structural eqv. wtd. DV","Num. interlocks", "Eig. Centrality", "Num. opp. coalitions","Employees","Revenue") 
N <- c(sum(summary(mod1b_lag)$df[1:2]),sum(summary(mod2b_lag)$df[1:2]),sum(summary(mod3b_lag)$df[1:2]),sum(summary(mod4b_lag)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(22,23); addtorow$command <- c(paste0(' '),paste0(''))
ptable(cbind(rownames(tab), tab), "output/lag_saturated.tex")

mods <- list(summary(mod1_fd)$coef[1:2,],summary(mod2_fd)$coef[1:2,],summary(mod3_fd)$coef[1:2,],summary(mod4_fd)$coef[1:2,])
vars <- c(rownames(mods[[1]])[2],rownames(mods[[2]])[2],rownames(mods[[3]])[2],rownames(mods[[4]])[2],"(Intercept)")
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9)] <- 
  c("CSOs, Interlock wtd. (fd)","CDP reporting, Interlock wtd. (fd)","Climate coalitions, Interlock wtd. (fd)","Climate lobbying, Interlock wtd. (fd)","Intercept") 
N <- c(sum(summary(mod1_fd)$df[1:2]),sum(summary(mod2_fd)$df[1:2]),sum(summary(mod3_fd)$df[1:2]),sum(summary(mod4_fd)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(10,11); addtorow$command <- c(paste0(' '),paste0(' '))
ptable(cbind(rownames(tab), tab), "output/fd_bivariate.tex")

mods <- list(summary(mod1a_fd)$coef[1:8,],summary(mod2a_fd)$coef[1:8,],summary(mod3a_fd)$coef[1:8,],summary(mod4a_fd)$coef[1:8,])
rownames(mods[[1]])[3] <- rownames(mods[[2]])[3] <- rownames(mods[[3]])[3] <- rownames(mods[[4]])[3] <- "Structural"
vars <- c(rownames(mods[[1]])[2],rownames(mods[[2]])[2],rownames(mods[[3]])[2],rownames(mods[[4]])[2:8],"(Intercept)")
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17,19,21)] <- 
  c("CSOs, Interlock wtd. (fd)","CDP reporting, Interlock wtd. (fd)","Climate coalitions, Interlock wtd. (fd)","Climate lobbying, Interlock wtd. (fd)",
    "Structural eqv. wtd. DV","Num. interlocks", "Eig. Centrality", "Num. opp. coalitions","Employees","Revenue","Intercept") 
N <- c(sum(summary(mod1a_fd)$df[1:2]),sum(summary(mod2a_fd)$df[1:2]),sum(summary(mod3a_fd)$df[1:2]),sum(summary(mod4a_fd)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(22,23); addtorow$command <- c(paste0(' '),paste0(' '))
ptable(cbind(rownames(tab), tab), "output/fd_covariates.tex")

# clean up
rems <- ls()[ls() %in% c("dat","cur","pdat","pcur","myround","ptable","dat_lag","cur_lag") == FALSE]
rm(list = rems)

# 5. Models with indirect interlocks -------------------------------------
mod1_indirect <- lm(cso_exists ~ cso_exists_indilwtd, data = cur)
mod1aa_indirect <- lm(cso_exists ~ cso_exists_indilwtd + cso_exists_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = cur)
mod1a_indirect <- lm(cso_exists ~ cso_exists_indilwtd + cso_exists_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = cur)
mod1b_indirect <- lm(cso_exists ~ cso_exists_indilwtd + cso_exists_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = pcur)

mod2_indirect <- lm(cdp_report ~ cdp_report_indilwtd, data = dat)
mod2aa_indirect <- lm(cdp_report ~ cdp_report_indilwtd + cdp_report_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat)
mod2a_indirect <- lm(cdp_report ~ cdp_report_indilwtd + cdp_report_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = dat)
mod2b_indirect <- lm(cdp_report ~ cdp_report_indilwtd + cdp_report_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = pdat)

mod3_indirect <- lm(I(numsupcoal>0) ~ numsupcoal_indilwtd + numsupcoal_sewtd, data = dat)
mod3aa_indirect <- lm(I(numsupcoal>0) ~ numsupcoal_indilwtd + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat)
mod3a_indirect <- lm(I(numsupcoal>0) ~ numsupcoal_indilwtd + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = dat)
mod3b_indirect <- lm(numsupcoal_d ~ numsupcoal_indilwtd + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = pdat)

mod4_indirect <- lm(lob ~ lob_indilwtd, data = dat)
mod4aa_indirect <- lm(lob ~ lob_indilwtd + lob_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat)
mod4a_indirect <- lm(lob ~ lob_indilwtd + lob_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = dat)
mod4b_indirect <- lm(lob ~ lob_indilwtd + lob_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = pdat)

mods <- list(summary(mod1_indirect)$coef[1:2,],summary(mod2_indirect)$coef[1:2,],summary(mod3_indirect)$coef[1:2,],summary(mod4_indirect)$coef[1:2,])
vars <- c(rownames(mods[[1]])[2],rownames(mods[[2]])[2],rownames(mods[[3]])[2],rownames(mods[[4]])[2],"(Intercept)")
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9)] <- 
  c("CSOs, Indirect interlock wtd.","CDP reporting, Indirect interlock wtd.","Climate coalitions, Indirect interlock wtd.","Climate lobbying, Indirect interlock wtd.","Intercept") 
N <- c(sum(summary(mod1_indirect)$df[1:2]),sum(summary(mod2_indirect)$df[1:2]),sum(summary(mod3_indirect)$df[1:2]),sum(summary(mod4_indirect)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(10,11); addtorow$command <- c(paste0(' '),paste0(' '))
ptable(cbind(rownames(tab), tab), "output/indirect_bivariate.tex")

mods <- list(summary(mod1aa_indirect)$coef[2:8,],summary(mod2aa_indirect)$coef[2:8,],summary(mod3aa_indirect)$coef[2:8,],summary(mod4aa_indirect)$coef[2:8,])
rownames(mods[[1]])[2] <- rownames(mods[[2]])[2] <- rownames(mods[[3]])[2] <- rownames(mods[[4]])[2] <- "Structural"
vars <- c(rownames(mods[[1]])[1],rownames(mods[[2]])[1],rownames(mods[[3]])[1],rownames(mods[[4]])[1],rownames(mods[[1]])[2:7])
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17,19)] <- 
  c("CSOs, Indirect interlock wtd.","CDP reporting, Indirect interlock wtd.","Climate coalitions, Indirect interlock wtd.","Climate lobbying, Indirect interlock wtd.",
    "Structural eqv. wtd. DV","Num. interlocks", "Eig. Centrality","Num. opp. coalitions","Employees","Revenue") 
N <- c(sum(summary(mod1aa_indirect)$df[1:2]),sum(summary(mod2aa_indirect)$df[1:2]),sum(summary(mod3aa_indirect)$df[1:2]),sum(summary(mod4aa_indirect)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(20,21); addtorow$command <- c(paste0(' '),paste0(''))
ptable(cbind(rownames(tab), tab), "output/indirect_indnotfe.tex")

mods <- list(summary(mod1a_indirect)$coef[2:8,],summary(mod2a_indirect)$coef[2:8,],summary(mod3a_indirect)$coef[2:8,],summary(mod4a_indirect)$coef[2:8,])
rownames(mods[[1]])[2] <- rownames(mods[[2]])[2] <- rownames(mods[[3]])[2] <- rownames(mods[[4]])[2] <- "Structural"
vars <- c(rownames(mods[[1]])[1],rownames(mods[[2]])[1],rownames(mods[[3]])[1],rownames(mods[[4]])[1],rownames(mods[[1]])[2:7])
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17,19)] <- 
  c("CSOs, Indirect interlock wtd.","CDP reporting, Indirect interlock wtd.","Climate coalitions, Indirect interlock wtd.","Climate lobbying, Indirect interlock wtd.",
    "Structural eqv. wtd. DV","Num. interlocks", "Eig. Centrality","Num. opp. coalitions","Employees","Revenue") 
N <- c(sum(summary(mod1a_indirect)$df[1:2]),sum(summary(mod2a_indirect)$df[1:2]),sum(summary(mod3a_indirect)$df[1:2]),sum(summary(mod4a_indirect)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(20,21); addtorow$command <- c(paste0(' '),paste0(''))
ptable(cbind(rownames(tab), tab), "output/indirect_indfe.tex")

mods <- list(summary(mod1b_indirect)$coef[2:8,],summary(mod2b_indirect)$coef[2:8,],summary(mod3b_indirect)$coef[2:8,],summary(mod4b_indirect)$coef[2:8,])
rownames(mods[[1]])[2] <- rownames(mods[[2]])[2] <- rownames(mods[[3]])[2] <- rownames(mods[[4]])[2] <- "Structural"
vars <- c(rownames(mods[[1]])[1],rownames(mods[[2]])[1],rownames(mods[[3]])[1],rownames(mods[[4]])[1],rownames(mods[[1]])[2:7])
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17,19)] <- 
  c("CSOs, Indirect interlock wtd.","CDP reporting, Indirect interlock wtd.","Climate coalitions, Indirect interlock wtd.","Climate lobbying, Indirect interlock wtd.",
    "Structural eqv. wtd. DV","Num. interlocks", "Eig. Centrality","Num. opp. coalitions","Employees","Revenue") 
N <- c(sum(summary(mod1b_indirect)$df[1:2]),sum(summary(mod2b_indirect)$df[1:2]),sum(summary(mod3b_indirect)$df[1:2]),sum(summary(mod4b_indirect)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(20,21); addtorow$command <- c(paste0(' '),paste0(''))
ptable(cbind(rownames(tab), tab), "output/indirect_saturated.tex")

# clean up
rems <- ls()[ls() %in% c("dat","cur","pdat","pcur","myround","ptable","dat_lag","cur_lag") == FALSE]
rm(list = rems)

# 6. Models with both indirect and direct interlocks -------------------------------------
mod1_indirect_directcontrol <- lm(cso_exists ~ cso_exists_indilwtd + cso_exists_ilwtd, data = cur)
mod1aa_indirect_directcontrol <- lm(cso_exists ~ cso_exists_indilwtd + cso_exists_ilwtd + cso_exists_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = cur)
mod1a_indirect_directcontrol <- lm(cso_exists ~ cso_exists_indilwtd + cso_exists_ilwtd + cso_exists_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = cur)
mod1b_indirect_directcontrol <- lm(cso_exists ~ cso_exists_indilwtd + cso_exists_ilwtd + cso_exists_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = pcur)

mod2_indirect_directcontrol <- lm(cdp_report ~ cdp_report_indilwtd + cdp_report_ilwtd, data = dat)
mod2aa_indirect_directcontrol <- lm(cdp_report ~ cdp_report_indilwtd + cdp_report_ilwtd + cdp_report_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat)
mod2a_indirect_directcontrol <- lm(cdp_report ~ cdp_report_indilwtd + cdp_report_ilwtd + cdp_report_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = dat)
mod2b_indirect_directcontrol <- lm(cdp_report ~ cdp_report_indilwtd + cdp_report_ilwtd + cdp_report_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = pdat)

mod3_indirect_directcontrol <- lm(I(numsupcoal>0) ~ numsupcoal_indilwtd + numsupcoal_ilwtd, data = dat)
mod3aa_indirect_directcontrol <- lm(I(numsupcoal>0) ~ numsupcoal_indilwtd + numsupcoal_ilwtd + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat)
mod3a_indirect_directcontrol <- lm(I(numsupcoal>0) ~ numsupcoal_indilwtd + numsupcoal_ilwtd + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = dat)
mod3b_indirect_directcontrol <- lm(numsupcoal_d ~ numsupcoal_indilwtd + numsupcoal_ilwtd + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = pdat)

mod4_indirect_directcontrol <- lm(lob ~ lob_indilwtd + lob_ilwtd, data = dat)
mod4aa_indirect_directcontrol <- lm(lob ~ lob_indilwtd + lob_ilwtd + lob_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat)
mod4a_indirect_directcontrol <- lm(lob ~ lob_indilwtd + lob_ilwtd + lob_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = dat)
mod4b_indirect_directcontrol <- lm(lob ~ lob_indilwtd + lob_ilwtd + lob_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = pdat)

mods <- list(summary(mod1_indirect_directcontrol)$coef[1:3,],summary(mod2_indirect_directcontrol)$coef[1:3,],summary(mod3_indirect_directcontrol)$coef[1:3,],summary(mod4_indirect_directcontrol)$coef[1:3,])
rownames(mods[[1]])[3] <- rownames(mods[[2]])[3] <- rownames(mods[[3]])[3] <- rownames(mods[[4]])[3] <- "Direct"
vars <- c(rownames(mods[[1]])[2],rownames(mods[[2]])[2],rownames(mods[[3]])[2],rownames(mods[[4]])[2:3],"(Intercept)")
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11)] <- 
  c("CSOs, Indirect interlock wtd.","CDP reporting, Indirect interlock wtd.","Climate coalitions, Indirect interlock wtd.","Climate lobbying, Indirect interlock wtd.","Direct interlock wtd. DV","Intercept") 
N <- c(sum(summary(mod1_indirect_directcontrol)$df[1:2]),sum(summary(mod2_indirect_directcontrol)$df[1:2]),sum(summary(mod3_indirect_directcontrol)$df[1:2]),sum(summary(mod4_indirect_directcontrol)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(12,13); addtorow$command <- c(paste0(' '),paste0(' '))
ptable(cbind(rownames(tab), tab), "output/indirect_directcontrol_bivariate.tex")

mods <- list(summary(mod1aa_indirect_directcontrol)$coef[2:9,],summary(mod2aa_indirect_directcontrol)$coef[2:9,],summary(mod3aa_indirect_directcontrol)$coef[2:9,],summary(mod4aa_indirect_directcontrol)$coef[2:9,])
rownames(mods[[1]])[2] <- rownames(mods[[2]])[2] <- rownames(mods[[3]])[2] <- rownames(mods[[4]])[2] <- "Direct"
rownames(mods[[1]])[3] <- rownames(mods[[2]])[3] <- rownames(mods[[3]])[3] <- rownames(mods[[4]])[3] <- "Structural"
vars <- c(rownames(mods[[1]])[1],rownames(mods[[2]])[1],rownames(mods[[3]])[1],rownames(mods[[4]])[1],rownames(mods[[1]])[2:8])
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17,19,21)] <- 
  c("CSOs, Indirect interlock wtd.","CDP reporting, Indirect interlock wtd.","Climate coalitions, Indirect interlock wtd.","Climate lobbying, Indirect interlock wtd.",
    "Direct interlock wtd. DV","Structural eqv. wtd. DV","Num. interlocks","Eig. Centrality","Num. opp. coalitions","Employees","Revenue") 
N <- c(sum(summary(mod1aa_indirect_directcontrol)$df[1:2]),sum(summary(mod2aa_indirect_directcontrol)$df[1:2]),sum(summary(mod3aa_indirect_directcontrol)$df[1:2]),sum(summary(mod4aa_indirect_directcontrol)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(22,23); addtorow$command <- c(paste0(' '),paste0(''))
ptable(cbind(rownames(tab), tab), "output/indirect_directcontrol_indnotfe.tex")

mods <- list(summary(mod1a_indirect_directcontrol)$coef[2:9,],summary(mod2a_indirect_directcontrol)$coef[2:9,],summary(mod3a_indirect_directcontrol)$coef[2:9,],summary(mod4a_indirect_directcontrol)$coef[2:9,])
rownames(mods[[1]])[2] <- rownames(mods[[2]])[2] <- rownames(mods[[3]])[2] <- rownames(mods[[4]])[2] <- "Direct"
rownames(mods[[1]])[3] <- rownames(mods[[2]])[3] <- rownames(mods[[3]])[3] <- rownames(mods[[4]])[3] <- "Structural"
vars <- c(rownames(mods[[1]])[1],rownames(mods[[2]])[1],rownames(mods[[3]])[1],rownames(mods[[4]])[1],rownames(mods[[1]])[2:8])
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17,19,21)] <- 
  c("CSOs, Indirect interlock wtd.","CDP reporting, Indirect interlock wtd.","Climate coalitions, Indirect interlock wtd.","Climate lobbying, Indirect interlock wtd.",
    "Direct interlock wtd. DV","Structural eqv. wtd. DV","Num. interlocks", "Eig. Centrality","Num. opp. coalitions","Employees","Revenue") 
N <- c(sum(summary(mod1a_indirect_directcontrol)$df[1:2]),sum(summary(mod2a_indirect_directcontrol)$df[1:2]),sum(summary(mod3a_indirect_directcontrol)$df[1:2]),sum(summary(mod4a_indirect_directcontrol)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(22,23); addtorow$command <- c(paste0(' '),paste0(''))
ptable(cbind(rownames(tab), tab), "output/indirect_directcontrol_indfe.tex")

mods <- list(summary(mod1b_indirect_directcontrol)$coef[2:9,],summary(mod2b_indirect_directcontrol)$coef[2:9,],summary(mod3b_indirect_directcontrol)$coef[2:9,],summary(mod4b_indirect_directcontrol)$coef[2:9,])
rownames(mods[[1]])[2] <- rownames(mods[[2]])[2] <- rownames(mods[[3]])[2] <- rownames(mods[[4]])[2] <- "Direct"
rownames(mods[[1]])[3] <- rownames(mods[[2]])[3] <- rownames(mods[[3]])[3] <- rownames(mods[[4]])[3] <- "Structural"
vars <- c(rownames(mods[[1]])[1],rownames(mods[[2]])[1],rownames(mods[[3]])[1],rownames(mods[[4]])[1],rownames(mods[[1]])[2:8])
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17,19,21)] <- 
  c("CSOs, Indirect interlock wtd.","CDP reporting, Indirect interlock wtd.","Climate coalitions, Indirect interlock wtd.","Climate lobbying, Indirect interlock wtd.",
    "Direct interlock wtd. DV","Structural eqv. wtd. DV","Num. interlocks","Eig. Centrality","Num. opp. coalitions","Employees","Revenue") 
N <- c(sum(summary(mod1b_indirect_directcontrol)$df[1:2]),sum(summary(mod2b_indirect_directcontrol)$df[1:2]),sum(summary(mod3b_indirect_directcontrol)$df[1:2]),sum(summary(mod4b_indirect_directcontrol)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(22,23); addtorow$command <- c(paste0(' '),paste0(''))
ptable(cbind(rownames(tab), tab), "output/indirect_directcontrol_saturated.tex")

# clean up
rems <- ls()[ls() %in% c("dat","cur","pdat","pcur","myround","ptable","dat_lag","cur_lag") == FALSE]
rm(list = rems)

# 7. Models with binary indirect interlocks -------------------------------------
mod1_indirect_alt <- lm(cso_exists ~ cso_exists_indilwtd_alt, data = cur)
mod1aa_indirect_alt <- lm(cso_exists ~ cso_exists_indilwtd_alt + cso_exists_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = cur)
mod1a_indirect_alt <- lm(cso_exists ~ cso_exists_indilwtd_alt + cso_exists_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = cur)
mod1b_indirect_alt <- lm(cso_exists ~ cso_exists_indilwtd_alt + cso_exists_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = pcur)

mod2_indirect_alt <- lm(cdp_report ~ cdp_report_indilwtd_alt, data = dat)
mod2aa_indirect_alt <- lm(cdp_report ~ cdp_report_indilwtd_alt + cdp_report_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat)
mod2a_indirect_alt <- lm(cdp_report ~ cdp_report_indilwtd_alt + cdp_report_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = dat)
mod2b_indirect_alt <- lm(cdp_report ~ cdp_report_indilwtd_alt + cdp_report_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = pdat)

mod3_indirect_alt <- lm(I(numsupcoal>0) ~ numsupcoal_indilwtd_alt + numsupcoal_sewtd, data = dat)
mod3aa_indirect_alt <- lm(I(numsupcoal>0) ~ numsupcoal_indilwtd_alt + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat)
mod3a_indirect_alt <- lm(I(numsupcoal>0) ~ numsupcoal_indilwtd_alt + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = dat)
mod3b_indirect_alt <- lm(numsupcoal_d ~ numsupcoal_indilwtd_alt + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = pdat)

mod4_indirect_alt <- lm(lob ~ lob_indilwtd_alt, data = dat)
mod4aa_indirect_alt <- lm(lob ~ lob_indilwtd_alt + lob_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat)
mod4a_indirect_alt <- lm(lob ~ lob_indilwtd_alt + lob_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = dat)
mod4b_indirect_alt <- lm(lob ~ lob_indilwtd_alt + lob_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = pdat)

mods <- list(summary(mod1_indirect_alt)$coef[1:2,],summary(mod2_indirect_alt)$coef[1:2,],summary(mod3_indirect_alt)$coef[1:2,],summary(mod4_indirect_alt)$coef[1:2,])
vars <- c(rownames(mods[[1]])[2],rownames(mods[[2]])[2],rownames(mods[[3]])[2],rownames(mods[[4]])[2],"(Intercept)")
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9)] <- 
  c("CSOs, Indirect interlock wtd. (binary, no direct ties)","CDP reporting, Indirect interlock wtd. (binary, no direct ties)","Climate coalitions, Indirect interlock wtd. (binary, no direct ties)","Climate lobbying, Indirect interlock wtd. (binary, no direct ties)","Intercept") 
N <- c(sum(summary(mod1_indirect_alt)$df[1:2]),sum(summary(mod2_indirect_alt)$df[1:2]),sum(summary(mod3_indirect_alt)$df[1:2]),sum(summary(mod4_indirect_alt)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(10,11); addtorow$command <- c(paste0(' '),paste0(' '))
ptable(cbind(rownames(tab), tab), "output/indirect_alt_bivariate.tex")

mods <- list(summary(mod1aa_indirect_alt)$coef[2:8,],summary(mod2aa_indirect_alt)$coef[2:8,],summary(mod3aa_indirect_alt)$coef[2:8,],summary(mod4aa_indirect_alt)$coef[2:8,])
rownames(mods[[1]])[2] <- rownames(mods[[2]])[2] <- rownames(mods[[3]])[2] <- rownames(mods[[4]])[2] <- "Structural"
vars <- c(rownames(mods[[1]])[1],rownames(mods[[2]])[1],rownames(mods[[3]])[1],rownames(mods[[4]])[1],rownames(mods[[1]])[2:7])
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17,19)] <- 
  c("CSOs, Indirect interlock wtd. (binary, no direct ties)","CDP reporting, Indirect interlock wtd. (binary, no direct ties)","Climate coalitions, Indirect interlock wtd. (binary, no direct ties)","Climate lobbying, Indirect interlock wtd. (binary, no direct ties)",
    "Structural eqv. wtd. DV","Num. interlocks", "Eig. Centrality","Num. opp. coalitions","Employees","Revenue") 
N <- c(sum(summary(mod1aa_indirect_alt)$df[1:2]),sum(summary(mod2aa_indirect_alt)$df[1:2]),sum(summary(mod3aa_indirect_alt)$df[1:2]),sum(summary(mod4aa_indirect_alt)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(20,21); addtorow$command <- c(paste0(' '),paste0(''))
ptable(cbind(rownames(tab), tab), "output/indirect_alt_indnotfe.tex")

mods <- list(summary(mod1a_indirect_alt)$coef[2:8,],summary(mod2a_indirect_alt)$coef[2:8,],summary(mod3a_indirect_alt)$coef[2:8,],summary(mod4a_indirect_alt)$coef[2:8,])
rownames(mods[[1]])[2] <- rownames(mods[[2]])[2] <- rownames(mods[[3]])[2] <- rownames(mods[[4]])[2] <- "Structural"
vars <- c(rownames(mods[[1]])[1],rownames(mods[[2]])[1],rownames(mods[[3]])[1],rownames(mods[[4]])[1],rownames(mods[[1]])[2:7])
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17,19)] <- 
  c("CSOs, Indirect interlock wtd. (binary, no direct ties)","CDP reporting, Indirect interlock wtd. (binary, no direct ties)","Climate coalitions, Indirect interlock wtd. (binary, no direct ties)","Climate lobbying, Indirect interlock wtd. (binary, no direct ties)",
    "Structural eqv. wtd. DV","Num. interlocks","Eig. Centrality","Num. opp. coalitions","Employees","Revenue") 
N <- c(sum(summary(mod1a_indirect_alt)$df[1:2]),sum(summary(mod2a_indirect_alt)$df[1:2]),sum(summary(mod3a_indirect_alt)$df[1:2]),sum(summary(mod4a_indirect_alt)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(20,21); addtorow$command <- c(paste0(' '),paste0(''))
ptable(cbind(rownames(tab), tab), "output/indirect_alt_indfe.tex")

mods <- list(summary(mod1b_indirect_alt)$coef[2:8,],summary(mod2b_indirect_alt)$coef[2:8,],summary(mod3b_indirect_alt)$coef[2:8,],summary(mod4b_indirect_alt)$coef[2:8,])
rownames(mods[[1]])[2] <- rownames(mods[[2]])[2] <- rownames(mods[[3]])[2] <- rownames(mods[[4]])[2] <- "Structural"
vars <- c(rownames(mods[[1]])[1],rownames(mods[[2]])[1],rownames(mods[[3]])[1],rownames(mods[[4]])[1],rownames(mods[[1]])[2:7])
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17,19)] <- 
  c("CSOs, Indirect interlock wtd. (binary, no direct ties)","CDP reporting, Indirect interlock wtd. (binary, no direct ties)","Climate coalitions, Indirect interlock wtd. (binary, no direct ties)","Climate lobbying, Indirect interlock wtd. (binary, no direct ties)",
    "Structural eqv. wtd. DV","Num. interlocks","Eig. Centrality","Num. opp. coalitions","Employees","Revenue") 
N <- c(sum(summary(mod1b_indirect_alt)$df[1:2]),sum(summary(mod2b_indirect_alt)$df[1:2]),sum(summary(mod3b_indirect_alt)$df[1:2]),sum(summary(mod4b_indirect_alt)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(20,21); addtorow$command <- c(paste0(' '),paste0(''))
ptable(cbind(rownames(tab), tab), "output/indirect_alt_saturated.tex")

# clean up
rems <- ls()[ls() %in% c("dat","cur","pdat","pcur","myround","ptable","dat_lag","cur_lag") == FALSE]
rm(list = rems)

# 8. Models with initial adoption only -------------------------------------
initial_dat_func <- function(dat, dv,to_pdata = "no"){
  colnames(dat)[which(colnames(dat)==dv)] <- "temp"
  
  dat <- dat %>%
    filter(is.na(temp)==F) %>%
    mutate(year = as.numeric(as.character(year))) %>%
    group_by(gvkey) %>%
    mutate(drop = if_else(year!=min(year) & lag(temp,1)>=1,1,0)) %>%
    mutate(drop = if_else(cumsum(drop)>=1,1,drop)) %>%
    filter(drop == 0) %>%
    mutate(year = as.factor(year))
  
  colnames(dat)[which(colnames(dat)=="temp")] <- dv
  
  if(to_pdata=="pdat"){
    # Create pdata.frame and demean by firm 
    pdat <- pdata.frame(dat,index = c("gvkey","year"))
    varlist <- colnames(pdat)[-which(colnames(pdat) %in% c("gvkey","conml","year","cik","cusip","naics","emissions_direct_intensity", "ind_year","ind"))]
    for(v in varlist){
      pdat[,v] <- Within(pdat[,v],effect = "individual",na.rm = T)
    }
    return(pdat)
  } else if (to_pdata == "pcur"){
    # Subset curtailed data frame
    cur <- dat[dat$year %in% 2004:2021,]
    # Create pdata.frame and demean by firm 
    pcur <- pdata.frame(cur,index = c("gvkey","year"))
    varlist <- colnames(pdat)[-which(colnames(pdat) %in% c("gvkey","conml","year","cik","cusip","naics","emissions_direct_intensity", "ind_year","ind"))]
    for(v in varlist){
      pcur[,v] <- Within(pcur[,v],effect = "individual",na.rm = T)
    }
    return(pcur)
  } else {
    return(dat)
  }
}

mod1_initial <- lm(cso_exists ~ cso_exists_ilwtd, data = initial_dat_func(cur, "cso_exists"))
mod1aa_initial <- lm(cso_exists ~ cso_exists_ilwtd + cso_exists_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = initial_dat_func(cur, "cso_exists"))
mod1a_initial <- lm(cso_exists ~ cso_exists_ilwtd + cso_exists_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = initial_dat_func(cur, "cso_exists"))
mod1b_initial <- lm(cso_exists ~ cso_exists_ilwtd + cso_exists_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = initial_dat_func(dat, "cso_exists",to_pdata = "pcur"))

mod2_initial <- lm(cdp_report ~ cdp_report_ilwtd, data = initial_dat_func(dat, "cdp_report"))
mod2aa_initial <- lm(cdp_report ~ cdp_report_ilwtd + cdp_report_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = initial_dat_func(dat, "cdp_report"))
mod2a_initial <- lm(cdp_report ~ cdp_report_ilwtd + cdp_report_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = initial_dat_func(dat, "cdp_report"))
mod2b_initial <- lm(cdp_report ~ cdp_report_ilwtd + cdp_report_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = initial_dat_func(dat, "cdp_report",to_pdata = "pdat"))

mod3_initial <- lm(I(numsupcoal>0) ~ numsupcoal_ilwtd, data = initial_dat_func(dat, "numsupcoal"))
mod3aa_initial <- lm(I(numsupcoal>0) ~ numsupcoal_ilwtd + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data =  initial_dat_func(dat, "numsupcoal"))
mod3a_initial <- lm(I(numsupcoal>0) ~ numsupcoal_ilwtd + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = initial_dat_func(dat, "numsupcoal"))
mod3b_initial <- lm(numsupcoal_d ~ numsupcoal_ilwtd + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = initial_dat_func(dat, "numsupcoal_d",to_pdata = "pdat"))

mod4_initial <- lm(lob ~ lob_ilwtd, data = initial_dat_func(dat, "lob"))
mod4aa_initial <- lm(lob ~ lob_ilwtd + lob_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = initial_dat_func(dat, "lob"))
mod4a_initial <- lm(lob ~ lob_ilwtd + lob_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = initial_dat_func(dat, "lob"))
mod4b_initial <- lm(lob ~ lob_ilwtd + lob_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = initial_dat_func(dat, "lob",to_pdata = "pdat"))

mods <- list(summary(mod1_initial)$coef[1:2,],summary(mod2_initial)$coef[1:2,],summary(mod3_initial)$coef[1:2,],summary(mod4_initial)$coef[1:2,])
vars <- c(rownames(mods[[1]])[2],rownames(mods[[2]])[2],rownames(mods[[3]])[2],rownames(mods[[4]])[2],"(Intercept)")
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9)] <- 
  c("CSOs, Interlock wtd.","CDP reporting, Interlock wtd.","Climate coalitions, Interlock wtd.","Climate lobbying, Interlock wtd.","Intercept") 
N <- c(sum(summary(mod1_initial)$df[1:2]),sum(summary(mod2_initial)$df[1:2]),sum(summary(mod3_initial)$df[1:2]),sum(summary(mod4_initial)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(10,11); addtorow$command <- c(paste0(' '),paste0(' '))
ptable(cbind(rownames(tab), tab), "output/initial_bivariate.tex")

mods <- list(summary(mod1aa_initial)$coef[2:8,],summary(mod2aa_initial)$coef[2:8,],summary(mod3aa_initial)$coef[2:8,],summary(mod4aa_initial)$coef[2:8,])
rownames(mods[[1]])[2] <- rownames(mods[[2]])[2] <- rownames(mods[[3]])[2] <- rownames(mods[[4]])[2] <- "Structural"
vars <- c(rownames(mods[[1]])[1],rownames(mods[[2]])[1],rownames(mods[[3]])[1],rownames(mods[[4]])[1],rownames(mods[[1]])[2:7])
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17,19)] <- 
  c("CSOs, Interlock wtd.","CDP reporting, Interlock wtd.","Climate coalitions, Interlock wtd.","Climate lobbying, Interlock wtd.",
    "Structural eqv. wtd. DV","Num. interlocks", "Eig. Centrality", "Num. opp. coalitions","Employees","Revenue") 
N <- c(sum(summary(mod1aa_initial)$df[1:2]),sum(summary(mod2aa_initial)$df[1:2]),sum(summary(mod3aa_initial)$df[1:2]),sum(summary(mod4aa_initial)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(20,21); addtorow$command <- c(paste0(' '),paste0(''))
ptable(cbind(rownames(tab), tab), "output/initial_indnotyearfe.tex")

mods <- list(summary(mod1a_initial)$coef[2:8,],summary(mod2a_initial)$coef[2:8,],summary(mod3a_initial)$coef[2:8,],summary(mod4a_initial)$coef[2:8,])
rownames(mods[[1]])[2] <- rownames(mods[[2]])[2] <- rownames(mods[[3]])[2] <- rownames(mods[[4]])[2] <- "Structural"
vars <- c(rownames(mods[[1]])[1],rownames(mods[[2]])[1],rownames(mods[[3]])[1],rownames(mods[[4]])[1],rownames(mods[[1]])[2:7])
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17,19)] <- 
  c("CSOs, Interlock wtd.","CDP reporting, Interlock wtd.","Climate coalitions, Interlock wtd.","Climate lobbying, Interlock wtd.",
    "Structural eqv. wtd. DV","Num. interlocks", "Eig. Centrality", "Num. opp. coalitions","Employees","Revenue") 
N <- c(sum(summary(mod1a_initial)$df[1:2]),sum(summary(mod2a_initial)$df[1:2]),sum(summary(mod3a_initial)$df[1:2]),sum(summary(mod4a_initial)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(20,21); addtorow$command <- c(paste0(' '),paste0(''))
ptable(cbind(rownames(tab), tab), "output/initial_indfe.tex")

mods <- list(summary(mod1b_initial)$coef[2:8,],summary(mod2b_initial)$coef[2:8,],summary(mod3b_initial)$coef[2:8,],summary(mod4b_initial)$coef[2:8,])
rownames(mods[[1]])[2] <- rownames(mods[[2]])[2] <- rownames(mods[[3]])[2] <- rownames(mods[[4]])[2] <- "Structural"
vars <- c(rownames(mods[[1]])[1],rownames(mods[[2]])[1],rownames(mods[[3]])[1],rownames(mods[[4]])[1],rownames(mods[[1]])[2:7])
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17,19)] <- 
  c("CSOs, Interlock wtd.","CDP reporting, Interlock wtd.","Climate coalitions, Interlock wtd.","Climate lobbying, Interlock wtd.",
    "Structural eqv. wtd. DV","Num. interlocks", "Eig. Centrality", "Num. opp. coalitions","Employees","Revenue") 
N <- c(sum(summary(mod1b_initial)$df[1:2]),sum(summary(mod2b_initial)$df[1:2]),sum(summary(mod3b_initial)$df[1:2]),sum(summary(mod4b_initial)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(20,21); addtorow$command <- c(paste0(' '),paste0(''))
ptable(cbind(rownames(tab), tab), "output/initial_saturated.tex")

# clean up
rems <- ls()[ls() %in% c("dat","cur","pdat","pcur","myround","ptable","dat_lag","cur_lag") == FALSE]
rm(list = rems)

# 9. Models interacting interlock-weighted variables with sector emissions intensity -------------------------------------
mod1_co2 <- lm(cso_exists ~ cso_exists_ilwtd*emissions_direct_intensity, data = cur)
mod1aa_co2 <- lm(cso_exists ~ cso_exists_ilwtd*emissions_direct_intensity + cso_exists_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = cur)
mod1a_co2 <- lm(cso_exists ~ cso_exists_ilwtd*emissions_direct_intensity + cso_exists_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = cur)

mod2_co2 <- lm(cdp_report ~ cdp_report_ilwtd*emissions_direct_intensity, data = dat)
mod2aa_co2 <- lm(cdp_report ~ cdp_report_ilwtd*emissions_direct_intensity + cdp_report_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat)
mod2a_co2 <- lm(cdp_report ~ cdp_report_ilwtd*emissions_direct_intensity + cdp_report_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = dat)

mod3_co2 <- lm(I(numsupcoal>0) ~ numsupcoal_ilwtd*emissions_direct_intensity, data = dat)
mod3aa_co2 <- lm(I(numsupcoal>0) ~ numsupcoal_ilwtd*emissions_direct_intensity + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data =  dat)
mod3a_co2 <- lm(I(numsupcoal>0) ~ numsupcoal_ilwtd*emissions_direct_intensity + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = dat)

mod4_co2 <- lm(lob ~ lob_ilwtd*emissions_direct_intensity, data = dat)
mod4aa_co2 <- lm(lob ~ lob_ilwtd*emissions_direct_intensity + lob_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat)
mod4a_co2 <- lm(lob ~ lob_ilwtd*emissions_direct_intensity + lob_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = dat)

# Calculate linear combinations
mod1_co2_lincom_est <- as.vector(attr(linearHypothesis(mod1_co2, c("cso_exists_ilwtd+cso_exists_ilwtd:emissions_direct_intensity"),),"value"))
mod1_co2_lincom_se <- sqrt(as.vector(attr(linearHypothesis(mod1_co2, c("cso_exists_ilwtd+cso_exists_ilwtd:emissions_direct_intensity")),"vcov")))
mod1_co2_lincom_pval <- linearHypothesis(mod1_co2, c("cso_exists_ilwtd+cso_exists_ilwtd:emissions_direct_intensity"))$`Pr(>F)`[2]
mod1aa_co2_lincom_est <- as.vector(attr(linearHypothesis(mod1aa_co2, c("cso_exists_ilwtd+cso_exists_ilwtd:emissions_direct_intensity")),"value"))
mod1aa_co2_lincom_se <- sqrt(as.vector(attr(linearHypothesis(mod1aa_co2, c("cso_exists_ilwtd+cso_exists_ilwtd:emissions_direct_intensity")),"vcov")))
mod1aa_co2_lincom_pval <- linearHypothesis(mod1aa_co2, c("cso_exists_ilwtd+cso_exists_ilwtd:emissions_direct_intensity"))$`Pr(>F)`[2]
mod1a_co2_lincom_est <- as.vector(attr(linearHypothesis(mod1a_co2, c("cso_exists_ilwtd+cso_exists_ilwtd:emissions_direct_intensity")),"value"))
mod1a_co2_lincom_se <- sqrt(as.vector(attr(linearHypothesis(mod1a_co2, c("cso_exists_ilwtd+cso_exists_ilwtd:emissions_direct_intensity")),"vcov")))
mod1a_co2_lincom_pval <- linearHypothesis(mod1a_co2, c("cso_exists_ilwtd+cso_exists_ilwtd:emissions_direct_intensity"))$`Pr(>F)`[2]

mod2_co2_lincom_est <- as.vector(attr(linearHypothesis(mod2_co2, c("cdp_report_ilwtd+cdp_report_ilwtd:emissions_direct_intensity")),"value"))
mod2_co2_lincom_se <- sqrt(as.vector(attr(linearHypothesis(mod2_co2, c("cdp_report_ilwtd+cdp_report_ilwtd:emissions_direct_intensity")),"vcov")))
mod2_co2_lincom_pval <- linearHypothesis(mod2_co2, c("cdp_report_ilwtd+cdp_report_ilwtd:emissions_direct_intensity"))$`Pr(>F)`[2]
mod2aa_co2_lincom_est <- as.vector(attr(linearHypothesis(mod2aa_co2, c("cdp_report_ilwtd+cdp_report_ilwtd:emissions_direct_intensity")),"value"))
mod2aa_co2_lincom_se <- sqrt(as.vector(attr(linearHypothesis(mod2aa_co2, c("cdp_report_ilwtd+cdp_report_ilwtd:emissions_direct_intensity")),"vcov")))
mod2aa_co2_lincom_pval <- linearHypothesis(mod2aa_co2, c("cdp_report_ilwtd+cdp_report_ilwtd:emissions_direct_intensity"))$`Pr(>F)`[2]
mod2a_co2_lincom_est <- as.vector(attr(linearHypothesis(mod2a_co2, c("cdp_report_ilwtd+cdp_report_ilwtd:emissions_direct_intensity")),"value"))
mod2a_co2_lincom_se <- sqrt(as.vector(attr(linearHypothesis(mod2a_co2, c("cdp_report_ilwtd+cdp_report_ilwtd:emissions_direct_intensity")),"vcov")))
mod2a_co2_lincom_pval <- linearHypothesis(mod2a_co2, c("cdp_report_ilwtd+cdp_report_ilwtd:emissions_direct_intensity"))$`Pr(>F)`[2]

mod3_co2_lincom_est <- as.vector(attr(linearHypothesis(mod3_co2, c("numsupcoal_ilwtd+numsupcoal_ilwtd:emissions_direct_intensity")),"value"))
mod3_co2_lincom_se <- sqrt(as.vector(attr(linearHypothesis(mod3_co2, c("numsupcoal_ilwtd+numsupcoal_ilwtd:emissions_direct_intensity")),"vcov")))
mod3_co2_lincom_pval <- linearHypothesis(mod3_co2, c("numsupcoal_ilwtd+numsupcoal_ilwtd:emissions_direct_intensity"))$`Pr(>F)`[2]
mod3aa_co2_lincom_est <- as.vector(attr(linearHypothesis(mod3aa_co2, c("numsupcoal_ilwtd+numsupcoal_ilwtd:emissions_direct_intensity")),"value"))
mod3aa_co2_lincom_se <- sqrt(as.vector(attr(linearHypothesis(mod3aa_co2, c("numsupcoal_ilwtd+numsupcoal_ilwtd:emissions_direct_intensity")),"vcov")))
mod3aa_co2_lincom_pval <- linearHypothesis(mod3aa_co2, c("numsupcoal_ilwtd+numsupcoal_ilwtd:emissions_direct_intensity"))$`Pr(>F)`[2]
mod3a_co2_lincom_est <- as.vector(attr(linearHypothesis(mod3a_co2, c("numsupcoal_ilwtd+numsupcoal_ilwtd:emissions_direct_intensity")),"value"))
mod3a_co2_lincom_se <- sqrt(as.vector(attr(linearHypothesis(mod3a_co2, c("numsupcoal_ilwtd+numsupcoal_ilwtd:emissions_direct_intensity")),"vcov")))
mod3a_co2_lincom_pval <- linearHypothesis(mod3a_co2, c("numsupcoal_ilwtd+numsupcoal_ilwtd:emissions_direct_intensity"))$`Pr(>F)`[2]

mod4_co2_lincom_est <- as.vector(attr(linearHypothesis(mod4_co2, c("lob_ilwtd+lob_ilwtd:emissions_direct_intensity")),"value"))
mod4_co2_lincom_se <- sqrt(as.vector(attr(linearHypothesis(mod4_co2, c("lob_ilwtd+lob_ilwtd:emissions_direct_intensity")),"vcov")))
mod4_co2_lincom_pval <- linearHypothesis(mod4_co2, c("lob_ilwtd+lob_ilwtd:emissions_direct_intensity"))$`Pr(>F)`[2]
mod4aa_co2_lincom_est <- as.vector(attr(linearHypothesis(mod4aa_co2, c("lob_ilwtd+lob_ilwtd:emissions_direct_intensity")),"value"))
mod4aa_co2_lincom_se <- sqrt(as.vector(attr(linearHypothesis(mod4aa_co2, c("lob_ilwtd+lob_ilwtd:emissions_direct_intensity")),"vcov")))
mod4aa_co2_lincom_pval <- linearHypothesis(mod4aa_co2, c("lob_ilwtd+lob_ilwtd:emissions_direct_intensity"))$`Pr(>F)`[2]
mod4a_co2_lincom_est <- as.vector(attr(linearHypothesis(mod4a_co2, c("lob_ilwtd+lob_ilwtd:emissions_direct_intensity")),"value"))
mod4a_co2_lincom_se <- sqrt(as.vector(attr(linearHypothesis(mod4a_co2, c("lob_ilwtd+lob_ilwtd:emissions_direct_intensity")),"vcov")))
mod4a_co2_lincom_pval <- linearHypothesis(mod4a_co2, c("lob_ilwtd+lob_ilwtd:emissions_direct_intensity"))$`Pr(>F)`[2]

mods <- list(summary(mod1_co2)$coef[1:4,],summary(mod2_co2)$coef[1:4,],summary(mod3_co2)$coef[1:4,],summary(mod4_co2)$coef[1:4,])
vars <- c(rownames(mods[[1]])[2],rownames(mods[[2]])[2],rownames(mods[[3]])[2],rownames(mods[[4]])[2],"emissions_direct_intensity",
          paste0(c(rownames(mods[[1]])[2],rownames(mods[[2]])[2],rownames(mods[[3]])[2],rownames(mods[[4]])[2]),":emissions_direct_intensity"),
          "(Intercept)")
# Populate regression table                                                                                                                                          
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
# Condense interaction rows into one row
interact_tab <- matrix(rep(NA,2*ncol(tab)),nrow = 2)
rownames(interact_tab) <- rep("Direct interlock wtd. DV $\\times$ CO$_2$ emissions intensity",2)
interact_idx <- seq(grep(rownames(tab),pattern = "[:]")[1],grep(rownames(tab),pattern = "[:]")[length(grep(rownames(tab),pattern = "[:]"))], by = 2) 
for(i in interact_idx){
  interact_tab[1:2,which(interact_idx==i)] <- tab[c(i,i+1),which(interact_idx==i)]
}
tab <- rbind(rbind(tab[1:(ncol(tab)*2+2),],interact_tab),tab[c((max(grep(rownames(tab),pattern = "[:]"))+1):nrow(tab)),])
# Fix row names
rownames(tab) <- c(rbind(vars[-grep(vars,pattern = "[:]")[-1]], "")); rownames(tab)[seq(1,nrow(tab),by=2)] <- 
  c("CSOs, Interlock wtd.","CDP reporting, Interlock wtd.","Climate coalitions, Interlock wtd.","Climate lobbying, Interlock wtd.","CO$_2$ emissions intensity","Direct interlock wtd. DV $\\times$ CO$_2$ emissions intensity","Intercept") 
# Add N
N <- c(sum(summary(mod1_co2)$df[1:2]),sum(summary(mod2_co2)$df[1:2]),sum(summary(mod3_co2)$df[1:2]),sum(summary(mod4_co2)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(10,11); addtorow$command <- c(paste0(' '),paste0(' '))
ptable(cbind(rownames(tab), tab), "output/co2_bivariate.tex")
# Save table with linear combinations
lc <- matrix(c(mod1_co2_lincom_est, mod2_co2_lincom_est, mod3_co2_lincom_est, mod4_co2_lincom_est,
               mod1_co2_lincom_se, mod2_co2_lincom_se, mod3_co2_lincom_se, mod4_co2_lincom_se),nrow = 2,byrow = T)
lc_pvals <- c(mod1_co2_lincom_pval,mod2_co2_lincom_pval,mod3_co2_lincom_pval,mod4_co2_lincom_pval)
stars <- rep("", ncol(lc));stars[lc_pvals<.1] <- "^{+}"; stars[lc_pvals < .05] <- "^{*}"; 
stars[lc_pvals < .01] <- "^{**}"; stars[lc_pvals < .001] <- "^{***}";
lc <- rbind(paste(myround(lc[1,]*100,2),stars, sep = ""), paste("(",myround(lc[2,]*100,2),")",sep = ""))
rownames(lc) <- c("Direct interlock wtd. DV effect: High CO$_2$ emissions intensity","")
addtorow <- list(); addtorow$pos <- list(1,2); addtorow$command <- c(paste0(' '),paste0(' '))
ptable(cbind(rownames(lc), lc), "output/co2_bivariate_lc.tex")

mods <- list(summary(mod1aa_co2)$coef[c(2:9,nrow(summary(mod1aa_co2)$coef)),],summary(mod2aa_co2)$coef[c(2:9,nrow(summary(mod2aa_co2)$coef)),],
             summary(mod3aa_co2)$coef[c(2:9,nrow(summary(mod3aa_co2)$coef)),],summary(mod4aa_co2)$coef[c(2:9,nrow(summary(mod4aa_co2)$coef)),])
rownames(mods[[1]])[3] <- rownames(mods[[2]])[3] <- rownames(mods[[3]])[3] <- rownames(mods[[4]])[3] <- "Structural"
vars <- c(rownames(mods[[1]])[1],rownames(mods[[2]])[1],rownames(mods[[3]])[1],rownames(mods[[4]])[1],"emissions_direct_intensity",
          paste0(c(rownames(mods[[1]])[1],rownames(mods[[2]])[1],rownames(mods[[3]])[1],rownames(mods[[4]])[1]),":emissions_direct_intensity"),
          rownames(mods[[1]])[3:8])
# Populate regression table                                                                                                                                          
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
# Condense interaction rows into one row
interact_tab <- matrix(rep(NA,2*ncol(tab)),nrow = 2)
rownames(interact_tab) <- rep("Direct interlock wtd. DV $\\times$ CO$_2$ emissions intensity",2)
interact_idx <- seq(grep(rownames(tab),pattern = "[:]")[1],grep(rownames(tab),pattern = "[:]")[length(grep(rownames(tab),pattern = "[:]"))], by = 2) 
for(i in interact_idx){
  interact_tab[1:2,which(interact_idx==i)] <- tab[c(i,i+1),which(interact_idx==i)]
}
tab <- rbind(rbind(tab[1:(ncol(tab)*2+2),],interact_tab),tab[c((max(grep(rownames(tab),pattern = "[:]"))+1):nrow(tab)),])
# Fix row names
rownames(tab) <- c(rbind(vars[-grep(vars,pattern = "[:]")[-1]], "")); rownames(tab)[seq(1,nrow(tab),by=2)] <- 
  c("CSOs, Interlock wtd.","CDP reporting, Interlock wtd.","Climate coalitions, Interlock wtd.","Climate lobbying, Interlock wtd.",
    "CO$_2$ emissions intensity","Direct interlock wtd. DV $\\times$ CO$_2$ emissions intensity",
    "Structural eqv. wtd. DV","Num. interlocks", "Eig. Centrality", "Num. opp. coalitions","Employees","Revenue") 
# Add N
N <- c(sum(summary(mod1aa_co2)$df[1:2]),sum(summary(mod2aa_co2)$df[1:2]),sum(summary(mod3aa_co2)$df[1:2]),sum(summary(mod4aa_co2)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(10,11); addtorow$command <- c(paste0(' '),paste0(' '))
ptable(cbind(rownames(tab), tab), "output/co2_indnotyearfe.tex")
# Save table with linear combinations
lc <- matrix(c(mod1aa_co2_lincom_est, mod2aa_co2_lincom_est, mod3aa_co2_lincom_est, mod4aa_co2_lincom_est,
               mod1aa_co2_lincom_se, mod2aa_co2_lincom_se, mod3aa_co2_lincom_se, mod4aa_co2_lincom_se),nrow = 2,byrow = T)
lc_pvals <- c(mod1aa_co2_lincom_pval,mod2aa_co2_lincom_pval,mod3aa_co2_lincom_pval,mod4aa_co2_lincom_pval)
stars <- rep("", ncol(lc));stars[lc_pvals<.1] <- "^{+}"; stars[lc_pvals < .05] <- "^{*}"; 
stars[lc_pvals < .01] <- "^{**}"; stars[lc_pvals < .001] <- "^{***}";
lc <- rbind(paste(myround(lc[1,]*100,2),stars, sep = ""), paste("(",myround(lc[2,]*100,2),")",sep = ""))
rownames(lc) <- c("Direct interlock wtd. DV effect: High CO$_2$ emissions intensity","")
addtorow <- list(); addtorow$pos <- list(1,2); addtorow$command <- c(paste0(' '),paste0(' '))
ptable(cbind(rownames(lc), lc), "output/co2_indnotyearfe_lc.tex")

mods <- list(summary(mod1a_co2)$coef[c(2:9,nrow(summary(mod1a_co2)$coef)),],summary(mod2a_co2)$coef[c(2:9,nrow(summary(mod2a_co2)$coef)),],
             summary(mod3a_co2)$coef[c(2:9,nrow(summary(mod3a_co2)$coef)),],summary(mod4a_co2)$coef[c(2:9,nrow(summary(mod4a_co2)$coef)),])
rownames(mods[[1]])[3] <- rownames(mods[[2]])[3] <- rownames(mods[[3]])[3] <- rownames(mods[[4]])[3] <- "Structural"
vars <- c(rownames(mods[[1]])[1],rownames(mods[[2]])[1],rownames(mods[[3]])[1],rownames(mods[[4]])[1],"emissions_direct_intensity",
          paste0(c(rownames(mods[[1]])[1],rownames(mods[[2]])[1],rownames(mods[[3]])[1],rownames(mods[[4]])[1]),":emissions_direct_intensity"),
          rownames(mods[[1]])[3:8])
# Populate regression table                                                                                                                                          
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
# Condense interaction rows into one row
interact_tab <- matrix(rep(NA,2*ncol(tab)),nrow = 2)
rownames(interact_tab) <- rep("Direct interlock wtd. DV $\\times$ CO$_2$ emissions intensity",2)
interact_idx <- seq(grep(rownames(tab),pattern = "[:]")[1],grep(rownames(tab),pattern = "[:]")[length(grep(rownames(tab),pattern = "[:]"))], by = 2) 
for(i in interact_idx){
  interact_tab[1:2,which(interact_idx==i)] <- tab[c(i,i+1),which(interact_idx==i)]
}
tab <- rbind(rbind(tab[1:(ncol(tab)*2+2),],interact_tab),tab[c((max(grep(rownames(tab),pattern = "[:]"))+1):nrow(tab)),])
# Fix row names
rownames(tab) <- c(rbind(vars[-grep(vars,pattern = "[:]")[-1]], "")); rownames(tab)[seq(1,nrow(tab),by=2)] <- 
  c("CSOs, Interlock wtd.","CDP reporting, Interlock wtd.","Climate coalitions, Interlock wtd.","Climate lobbying, Interlock wtd.",
    "CO$_2$ emissions intensity","Direct interlock wtd. DV $\\times$ CO$_2$ emissions intensity",
    "Structural eqv. wtd. DV","Num. interlocks", "Eig. Centrality", "Num. opp. coalitions","Employees","Revenue") 
# Add N
N <- c(sum(summary(mod1a_co2)$df[1:2]),sum(summary(mod2a_co2)$df[1:2]),sum(summary(mod3a_co2)$df[1:2]),sum(summary(mod4a_co2)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(10,11); addtorow$command <- c(paste0(' '),paste0(' '))
ptable(cbind(rownames(tab), tab), "output/co2_indfe.tex")
# Save table with linear combinations
lc <- matrix(c(mod1a_co2_lincom_est, mod2a_co2_lincom_est, mod3a_co2_lincom_est, mod4a_co2_lincom_est,
               mod1a_co2_lincom_se, mod2a_co2_lincom_se, mod3a_co2_lincom_se, mod4a_co2_lincom_se),nrow = 2,byrow = T)
lc_pvals <- c(mod1a_co2_lincom_pval,mod2a_co2_lincom_pval,mod3a_co2_lincom_pval,mod4a_co2_lincom_pval)
stars <- rep("", ncol(lc));stars[lc_pvals<.1] <- "^{+}"; stars[lc_pvals < .05] <- "^{*}"; 
stars[lc_pvals < .01] <- "^{**}"; stars[lc_pvals < .001] <- "^{***}";
lc <- rbind(paste(myround(lc[1,]*100,2),stars, sep = ""), paste("(",myround(lc[2,]*100,2),")",sep = ""))
rownames(lc) <- c("Direct interlock wtd. DV effect: High CO$_2$ emissions intensity","")
addtorow <- list(); addtorow$pos <- list(1,2); addtorow$command <- c(paste0(' '),paste0(' '))
ptable(cbind(rownames(lc), lc), "output/co2_indfe_lc.tex")

# clean up
rems <- ls()[ls() %in% c("dat","cur","pdat","pcur","myround","ptable","dat_lag","cur_lag") == FALSE]
rm(list = rems)

# 10. Models with pro-climate lobbying firm interlocks -------------------------------------
mod4_lobpro <- lm(lob ~ lobsup_ilwtd, data = dat)
mod4aa_lobpro <- lm(lob ~ lobsup_ilwtd + lobsup_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat)
mod4a_lobpro <- lm(lob ~ lobsup_ilwtd + lobsup_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = dat)
mod4b_lobpro <- lm(lob ~ lobsup_ilwtd + lobsup_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = pdat)
mod4_glm_lobpro <- glm(lob ~ lobsup_ilwtd, data = dat, family = "binomial")
mod4aa_glm_lobpro <- glm(lob ~ lobsup_ilwtd + lobsup_sewtd + num_interlocks + evc + numoppcoal + emp + revt, data = dat, family = "binomial")

mod4_lobopp <- lm(I(lob & numoppcoal >0) ~ lobopp_ilwtd, data = dat)
mod4aa_lobopp <- lm(I(lob & numoppcoal >0) ~ lobopp_ilwtd + lobopp_sewtd + num_interlocks + evc + numsupcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat)
mod4a_lobopp <- lm(I(lob & numoppcoal >0) ~ lobopp_ilwtd + lobopp_sewtd + num_interlocks + evc + numsupcoal + emp + revt + as.factor(ind_year), data = dat)
mod4b_lobopp <- lm(I(lob & numoppcoal >0) ~ lobopp_ilwtd + lobopp_sewtd + num_interlocks + evc + numsupcoal + emp + revt + as.factor(ind_year), data = pdat)
mod4_glm_lobopp <- glm(I(lob & numoppcoal >0) ~ lobopp_ilwtd, data = dat, family = "binomial")
mod4aa_glm_lobopp <- glm(I(lob & numoppcoal >0) ~ lobopp_ilwtd + lobopp_sewtd + num_interlocks + evc + numsupcoal + emp + revt, data = dat, family = "binomial")

mod4_lobpro <- lm(I(lob & numsupcoal >0) ~ lob_ilwtd, data = dat)
mod4aa_lobpro <- lm(I(lob & numsupcoal >0) ~ lobsup_ilwtd + lobsup_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat)
mod4a_lobpro <- lm(I(lob & numsupcoal >0) ~ lobsup_ilwtd + lobsup_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = dat)
mod4b_lobpro <- lm(I(lob & numsupcoal >0) ~ lobsup_ilwtd + lobsup_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = pdat)
mod4_glm_lobpro <- glm(I(lob & numsupcoal >0) ~ lobsup_ilwtd, data = dat, family = "binomial")
mod4aa_glm_lobpro <- glm(I(lob & numsupcoal >0) ~ lobsup_ilwtd + lobsup_sewtd + num_interlocks + evc + numoppcoal + emp + revt, data = dat, family = "binomial")

mod4_lobopp_all <- lm(lob ~ lobopp_ilwtd, data = dat)
mod4aa_lobopp_all <- lm(lob ~ lobopp_ilwtd + lobopp_sewtd + num_interlocks + evc + numsupcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat)
mod4a_lobopp_all <- lm(lob ~ lobopp_ilwtd + lobopp_sewtd + num_interlocks + evc + numsupcoal + emp + revt + as.factor(ind_year), data = dat)
mod4b_lobopp_all <- lm(lob ~ lobopp_ilwtd + lobopp_sewtd + num_interlocks + evc + numsupcoal + emp + revt + as.factor(ind_year), data = pdat)
mod4_glm_lobopp_all <- glm(lob ~ lobopp_ilwtd, data = dat, family = "binomial")
mod4aa_glm_lobopp_all <- glm(lob ~ lobopp_ilwtd + lobopp_sewtd + num_interlocks + evc + numsupcoal + emp + revt, data = dat, family = "binomial")

mod4_lobpro_all <- lm(lob ~ lob_ilwtd, data = dat)
mod4aa_lobpro_all <- lm(lob ~ lobsup_ilwtd + lobsup_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat)
mod4a_lobpro_all <- lm(lob ~ lobsup_ilwtd + lobsup_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = dat)
mod4b_lobpro_all <- lm(lob ~ lobsup_ilwtd + lobsup_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = pdat)
mod4_glm_lobpro_all <- glm(lob ~ lobsup_ilwtd, data = dat, family = "binomial")
mod4aa_glm_lobpro_all <- glm(lob ~ lobsup_ilwtd + lobsup_sewtd + num_interlocks + evc + numoppcoal + emp + revt, data = dat, family = "binomial")

coefs_LOBPRO <- c(summary(mod4_lobpro)$coef[2,1],summary(mod4aa_lobpro)$coef[2,1],summary(mod4a_lobpro)$coef[2,1],summary(mod4b_lobpro)$coef[2,1],summary(mod4_glm_lobpro)$coef[2,1],summary(mod4aa_glm_lobpro)$coef[2,1])
se_LOBPRO <- c(summary(mod4_lobpro)$coef[2,2],summary(mod4aa_lobpro)$coef[2,2],summary(mod4a_lobpro)$coef[2,2],summary(mod4b_lobpro)$coef[2,2],summary(mod4_glm_lobpro)$coef[2,2],summary(mod4aa_glm_lobpro)$coef[2,2])
coefs_LOBPRO[1:4] <- 100*coefs_LOBPRO[1:4]; se_LOBPRO[1:4] <- 100*se_LOBPRO[1:4]
se_LOBPRO <- paste("(", myround(se_LOBPRO,2), ")", sep = "")
p_LOBPRO <- c(summary(mod4_lobpro)$coef[2,4],summary(mod4aa_lobpro)$coef[2,4],summary(mod4a_lobpro)$coef[2,4],summary(mod4b_lobpro)$coef[2,4],summary(mod4_glm_lobpro)$coef[2,4],summary(mod4aa_glm_lobpro)$coef[2,4])
stars <- rep("", length(coefs_LOBPRO)); stars[p_LOBPRO < .1] <- "^{+}"; stars[p_LOBPRO < .05] <- "^{*}"; stars[p_LOBPRO < .01] <- "^{**}"; stars[p_LOBPRO < .001] <- "^{***}";  
coefs_LOBPRO <- paste(myround(coefs_LOBPRO,2), stars, sep = "")

addtorow <- list(); addtorow$pos <- list(1); addtorow$command <- c("")
ptable(rbind(c("Climate lobbying (Pro), Interlock wtd.", coefs_LOBPRO),c("", se_LOBPRO)), "output/tab_LOBPRO.tex")

coefs_LOBPRO_all <- c(summary(mod4_lobpro_all)$coef[2,1],summary(mod4aa_lobpro_all)$coef[2,1],summary(mod4a_lobpro_all)$coef[2,1],summary(mod4b_lobpro_all)$coef[2,1],summary(mod4_glm_lobpro_all)$coef[2,1],summary(mod4aa_glm_lobpro_all)$coef[2,1])
se_LOBPRO_all <- c(summary(mod4_lobpro_all)$coef[2,2],summary(mod4aa_lobpro_all)$coef[2,2],summary(mod4a_lobpro_all)$coef[2,2],summary(mod4b_lobpro_all)$coef[2,2],summary(mod4_glm_lobpro_all)$coef[2,2],summary(mod4aa_glm_lobpro_all)$coef[2,2])
coefs_LOBPRO_all[1:4] <- 100*coefs_LOBPRO_all[1:4]; se_LOBPRO_all[1:4] <- 100*se_LOBPRO_all[1:4]
se_LOBPRO_all <- paste("(", myround(se_LOBPRO_all,2), ")", sep = "")
p_LOBPRO_all <- c(summary(mod4_lobpro_all)$coef[2,4],summary(mod4aa_lobpro_all)$coef[2,4],summary(mod4a_lobpro_all)$coef[2,4],summary(mod4b_lobpro_all)$coef[2,4],summary(mod4_glm_lobpro_all)$coef[2,4],summary(mod4aa_glm_lobpro_all)$coef[2,4])
stars <- rep("", length(coefs_LOBPRO_all)); stars[p_LOBPRO_all < .1] <- "^{+}"; stars[p_LOBPRO_all < .05] <- "^{*}"; stars[p_LOBPRO_all < .01] <- "^{**}"; stars[p_LOBPRO_all < .001] <- "^{***}";  
coefs_LOBPRO_all <- paste(myround(coefs_LOBPRO_all,2), stars, sep = "")

addtorow <- list(); addtorow$pos <- list(1); addtorow$command <- c("")
ptable(rbind(c("Climate lobbying (Pro), Interlock wtd.", coefs_LOBPRO_all),c("", se_LOBPRO_all)), "output/tab_LOBPRO_all.tex")

coefs_LOBOPP <- c(summary(mod4_lobopp)$coef[2,1],summary(mod4aa_lobopp)$coef[2,1],summary(mod4a_lobopp)$coef[2,1],summary(mod4b_lobopp)$coef[2,1],summary(mod4_glm_lobopp)$coef[2,1],summary(mod4aa_glm_lobopp)$coef[2,1])
se_LOBOPP <- c(summary(mod4_lobopp)$coef[2,2],summary(mod4aa_lobopp)$coef[2,2],summary(mod4a_lobopp)$coef[2,2],summary(mod4b_lobopp)$coef[2,2],summary(mod4_glm_lobopp)$coef[2,2],summary(mod4aa_glm_lobopp)$coef[2,2])
coefs_LOBOPP[1:4] <- 100*coefs_LOBOPP[1:4]; se_LOBOPP[1:4] <- 100*se_LOBOPP[1:4]
se_LOBOPP <- paste("(", myround(se_LOBOPP,2), ")", sep = "")
p_LOBOPP <- c(summary(mod4_lobopp)$coef[2,4],summary(mod4aa_lobopp)$coef[2,4],summary(mod4a_lobopp)$coef[2,4],summary(mod4b_lobopp)$coef[2,4],summary(mod4_glm_lobopp)$coef[2,4],summary(mod4aa_glm_lobopp)$coef[2,4])
stars <- rep("", length(coefs_LOBOPP)); stars[p_LOBOPP < .1] <- "^{+}"; stars[p_LOBOPP < .05] <- "^{*}"; stars[p_LOBOPP < .01] <- "^{**}"; stars[p_LOBOPP < .001] <- "^{***}";  
coefs_LOBOPP <- paste(myround(coefs_LOBOPP,2), stars, sep = "")
# coefs_LOBOPP[5:6] <- "\\multicolumn{1}{c}{-}"; se_LOBOPP[5:6] <- "\\multicolumn{1}{c}{-}"

addtorow <- list(); addtorow$pos <- list(1); addtorow$command <- c("")
ptable(rbind(c("Climate lobbying (Opp), Interlock wtd.", coefs_LOBOPP),c("", se_LOBOPP)), "output/tab_LOBOPP.tex")

coefs_LOBOPP_all <- c(summary(mod4_lobopp_all)$coef[2,1],summary(mod4aa_lobopp_all)$coef[2,1],summary(mod4a_lobopp_all)$coef[2,1],summary(mod4b_lobopp_all)$coef[2,1],summary(mod4_glm_lobopp_all)$coef[2,1],summary(mod4aa_glm_lobopp_all)$coef[2,1])
se_LOBOPP_all <- c(summary(mod4_lobopp_all)$coef[2,2],summary(mod4aa_lobopp_all)$coef[2,2],summary(mod4a_lobopp_all)$coef[2,2],summary(mod4b_lobopp_all)$coef[2,2],summary(mod4_glm_lobopp_all)$coef[2,2],summary(mod4aa_glm_lobopp_all)$coef[2,2])
coefs_LOBOPP_all[1:4] <- 100*coefs_LOBOPP_all[1:4]; se_LOBOPP_all[1:4] <- 100*se_LOBOPP_all[1:4]
se_LOBOPP_all <- paste("(", myround(se_LOBOPP_all,2), ")", sep = "")
p_LOBOPP_all <- c(summary(mod4_lobopp_all)$coef[2,4],summary(mod4aa_lobopp_all)$coef[2,4],summary(mod4a_lobopp_all)$coef[2,4],summary(mod4b_lobopp_all)$coef[2,4],summary(mod4_glm_lobopp_all)$coef[2,4],summary(mod4aa_glm_lobopp_all)$coef[2,4])
stars <- rep("", length(coefs_LOBOPP_all)); stars[p_LOBOPP_all < .1] <- "^{+}"; stars[p_LOBOPP_all < .05] <- "^{*}"; stars[p_LOBOPP_all < .01] <- "^{**}"; stars[p_LOBOPP_all < .001] <- "^{***}";  
coefs_LOBOPP_all <- paste(myround(coefs_LOBOPP_all,2), stars, sep = "")

addtorow <- list(); addtorow$pos <- list(1); addtorow$command <- c("")
ptable(rbind(c("Climate lobbying (Opp), Interlock wtd.", coefs_LOBOPP_all),c("", se_LOBOPP_all)), "output/tab_LOBOPP_all.tex")

# clean up
rems <- ls()[ls() %in% c("dat","cur","pdat","pcur","myround","ptable","dat_lag","cur_lag") == FALSE]
rm(list = rems)

# 11. Models interacting interlock-weighted variables with time -------------------------------------
mod1_time <- lm(cso_exists ~ cso_exists_ilwtd*I(year-2004), data = cur)
mod1aa_time <- lm(cso_exists ~ cso_exists_ilwtd + cso_exists_ilwtd:I(as.numeric(as.character(year))-2004) + cso_exists_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = cur)
mod1a_time <- lm(cso_exists ~ cso_exists_ilwtd + cso_exists_ilwtd:I(as.numeric(as.character(year))-2004) + cso_exists_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = cur)
mod1b_time <- lm(cso_exists ~ cso_exists_ilwtd + cso_exists_ilwtd:I(as.numeric(as.character(year))-2004) + cso_exists_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = pcur)

mod2_time <- lm(cdp_report ~ cdp_report_ilwtd*I(year-2007), data = dat)
mod2aa_time <- lm(cdp_report ~ cdp_report_ilwtd + cdp_report_ilwtd:I(as.numeric(as.character(year))-2007) + cdp_report_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat)
mod2a_time <- lm(cdp_report ~ cdp_report_ilwtd + cdp_report_ilwtd:I(as.numeric(as.character(year))-2007) + cdp_report_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = dat)
mod2b_time <- lm(cdp_report ~ cdp_report_ilwtd + cdp_report_ilwtd:I(as.numeric(as.character(year))-2007) + cdp_report_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = pdat)

mod3_time <- lm(I(numsupcoal>0) ~ numsupcoal_ilwtd*I(year-1996), data = dat)
mod3aa_time <- lm(I(numsupcoal>0) ~ numsupcoal_ilwtd + numsupcoal_ilwtd:I(as.numeric(as.character(year))-1996) + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat)
mod3a_time <- lm(I(numsupcoal>0) ~ numsupcoal_ilwtd + numsupcoal_ilwtd:I(as.numeric(as.character(year))-1996) + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = dat)
mod3b_time <- lm(numsupcoal_d ~ numsupcoal_ilwtd + numsupcoal_ilwtd:I(as.numeric(as.character(year))-1996) + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = pdat)

mod4_time <- lm(lob ~ lob_ilwtd*I(year-2001), data = dat)
mod4aa_time <- lm(lob ~ lob_ilwtd + lob_ilwtd:I(as.numeric(as.character(year))-2001) + lob_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind) + as.factor(year), data = dat)
mod4a_time <- lm(lob ~ lob_ilwtd + lob_ilwtd:I(as.numeric(as.character(year))-2001) + lob_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = dat)
mod4b_time <- lm(lob ~ lob_ilwtd + lob_ilwtd:I(as.numeric(as.character(year))-2001) + lob_sewtd + num_interlocks + evc + numoppcoal + emp + revt + as.factor(ind_year), data = pdat)

mods <- list(summary(mod1_time)$coef[1:4,],summary(mod2_time)$coef[1:4,],summary(mod3_time)$coef[1:4,],summary(mod4_time)$coef[1:4,])
rownames(mods[[1]])[3] <- rownames(mods[[2]])[3] <- rownames(mods[[3]])[3] <- rownames(mods[[4]])[3] <- "Years after base year"
vars <- c(rownames(mods[[1]])[2],rownames(mods[[2]])[2],rownames(mods[[3]])[2],rownames(mods[[4]])[2],"Years after base year",
          rownames(mods[[1]])[4],rownames(mods[[2]])[4],rownames(mods[[3]])[4],rownames(mods[[4]])[4],"(Intercept)")
# Populate regression table                                                                                                                                          
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
# Condense interaction rows into one row
interact_tab <- matrix(rep(NA,2*ncol(tab)),nrow = 2)
rownames(interact_tab) <- rep("Direct interlock wtd. DV $\\times$ Years after base year",2)
interact_idx <- seq(grep(rownames(tab),pattern = "[:]")[1],grep(rownames(tab),pattern = "[:]")[length(grep(rownames(tab),pattern = "[:]"))], by = 2) 
for(i in interact_idx){
  interact_tab[1:2,which(interact_idx==i)] <- tab[c(i,i+1),which(interact_idx==i)]
}
tab <- rbind(rbind(tab[1:(ncol(tab)*2+2),],interact_tab),tab[c((max(grep(rownames(tab),pattern = "[:]"))+1):nrow(tab)),])
# Fix row names
rownames(tab) <- c(rbind(vars[-grep(vars,pattern = "[:]")[-1]], "")); rownames(tab)[seq(1,nrow(tab),by=2)] <- 
  c("CSOs, Interlock wtd.","CDP reporting, Interlock wtd.","Climate coalitions, Interlock wtd.","Climate lobbying, Interlock wtd.","CO$_2$ emissions intensity","Direct interlock wtd. DV $\\times$ Years after base year","Intercept") 
# Add N
N <- c(sum(summary(mod1_time)$df[1:2]),sum(summary(mod2_time)$df[1:2]),sum(summary(mod3_time)$df[1:2]),sum(summary(mod4_time)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(10,11); addtorow$command <- c(paste0(' '),paste0(' '))
ptable(cbind(rownames(tab), tab), "output/time_bivariate.tex")

mods <- list(summary(mod1aa_time)$coef[c(2:8,nrow(summary(mod1aa_time)$coef)),],summary(mod2aa_time)$coef[c(2:8,nrow(summary(mod2aa_time)$coef)),],
             summary(mod3aa_time)$coef[c(2:8,nrow(summary(mod3aa_time)$coef)),],summary(mod4aa_time)$coef[c(2:8,nrow(summary(mod4aa_time)$coef)),])
rownames(mods[[1]])[2] <- rownames(mods[[2]])[2] <- rownames(mods[[3]])[2] <- rownames(mods[[4]])[2] <- "Structural"
vars <- c(rownames(mods[[1]])[1],rownames(mods[[2]])[1],rownames(mods[[3]])[1],rownames(mods[[4]])[1],
          rownames(mods[[1]])[8],rownames(mods[[2]])[8],rownames(mods[[3]])[8],rownames(mods[[4]])[8], rownames(mods[[1]])[2:7])
# Populate regression table                                                                                                                                          
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
# Condense interaction rows into one row
interact_tab <- matrix(rep(NA,2*ncol(tab)),nrow = 2)
rownames(interact_tab) <- rep("Direct interlock wtd. DV $\\times$ Years after base year",2)
interact_idx <- seq(grep(rownames(tab),pattern = "[:]")[1],grep(rownames(tab),pattern = "[:]")[length(grep(rownames(tab),pattern = "[:]"))], by = 2) 
for(i in interact_idx){
  interact_tab[1:2,which(interact_idx==i)] <- tab[c(i,i+1),which(interact_idx==i)]
}
tab <- rbind(rbind(tab[1:(ncol(tab)*2),],interact_tab),tab[c((max(grep(rownames(tab),pattern = "[:]"))+1):nrow(tab)),])
# Fix row names
rownames(tab) <- c(rbind(vars[-grep(vars,pattern = "[:]")[-1]], "")); 
rownames(tab)[seq(1,nrow(tab),by=2)] <- 
  c("CSOs, Interlock wtd.","CDP reporting, Interlock wtd.","Climate coalitions, Interlock wtd.","Climate lobbying, Interlock wtd.",
    "Direct interlock wtd. DV $\\times$ Years after base year",
    "Structural eqv. wtd. DV","Num. interlocks", "Eig. Centrality", "Num. opp. coalitions","Employees","Revenue") 
# Add N
N <- c(sum(summary(mod1aa_time)$df[1:2]),sum(summary(mod2aa_time)$df[1:2]),sum(summary(mod3aa_time)$df[1:2]),sum(summary(mod4aa_time)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(10,11); addtorow$command <- c(paste0(' '),paste0(' '))
ptable(cbind(rownames(tab), tab), "output/time_indnotyearfe.tex")

mods <- list(summary(mod1a_time)$coef[c(2:8,nrow(summary(mod1a_time)$coef)),],summary(mod2a_time)$coef[c(2:8,nrow(summary(mod2a_time)$coef)),],
             summary(mod3a_time)$coef[c(2:8,nrow(summary(mod3a_time)$coef)),],summary(mod4a_time)$coef[c(2:8,nrow(summary(mod4a_time)$coef)),])
rownames(mods[[1]])[2] <- rownames(mods[[2]])[2] <- rownames(mods[[3]])[2] <- rownames(mods[[4]])[2] <- "Structural"
vars <- c(rownames(mods[[1]])[1],rownames(mods[[2]])[1],rownames(mods[[3]])[1],rownames(mods[[4]])[1],
          rownames(mods[[1]])[8],rownames(mods[[2]])[8],rownames(mods[[3]])[8],rownames(mods[[4]])[8], rownames(mods[[1]])[2:7])
# Populate regression table                                                                                                                                          
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
# Condense interaction rows into one row
interact_tab <- matrix(rep(NA,2*ncol(tab)),nrow = 2)
rownames(interact_tab) <- rep("Direct interlock wtd. DV $\\times$ Years after base year",2)
interact_idx <- seq(grep(rownames(tab),pattern = "[:]")[1],grep(rownames(tab),pattern = "[:]")[length(grep(rownames(tab),pattern = "[:]"))], by = 2) 
for(i in interact_idx){
  interact_tab[1:2,which(interact_idx==i)] <- tab[c(i,i+1),which(interact_idx==i)]
}
tab <- rbind(rbind(tab[1:(ncol(tab)*2),],interact_tab),tab[c((max(grep(rownames(tab),pattern = "[:]"))+1):nrow(tab)),])
# Fix row names
rownames(tab) <- c(rbind(vars[-grep(vars,pattern = "[:]")[-1]], "")); 
rownames(tab)[seq(1,nrow(tab),by=2)] <- 
  c("CSOs, Interlock wtd.","CDP reporting, Interlock wtd.","Climate coalitions, Interlock wtd.","Climate lobbying, Interlock wtd.",
    "Direct interlock wtd. DV $\\times$ Years after base year",
    "Structural eqv. wtd. DV","Num. interlocks", "Eig. Centrality", "Num. opp. coalitions","Employees","Revenue") 
# Add N
N <- c(sum(summary(mod1a_time)$df[1:2]),sum(summary(mod2a_time)$df[1:2]),sum(summary(mod3a_time)$df[1:2]),sum(summary(mod4a_time)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(10,11); addtorow$command <- c(paste0(' '),paste0(' '))
ptable(cbind(rownames(tab), tab), "output/time_indyearfe.tex")

mods <- list(summary(mod1b_time)$coef[c(2:8,nrow(summary(mod1b_time)$coef)),],summary(mod2b_time)$coef[c(2:8,nrow(summary(mod2b_time)$coef)),],
             summary(mod3b_time)$coef[c(2:8,nrow(summary(mod3b_time)$coef)),],summary(mod4b_time)$coef[c(2:8,nrow(summary(mod4b_time)$coef)),])
rownames(mods[[1]])[2] <- rownames(mods[[2]])[2] <- rownames(mods[[3]])[2] <- rownames(mods[[4]])[2] <- "Structural"
vars <- c(rownames(mods[[1]])[1],rownames(mods[[2]])[1],rownames(mods[[3]])[1],rownames(mods[[4]])[1],
          rownames(mods[[1]])[8],rownames(mods[[2]])[8],rownames(mods[[3]])[8],rownames(mods[[4]])[8], rownames(mods[[1]])[2:7])
# Populate regression table                                                                                                                                          
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1]*100,2), stars, sep = ""), paste("(", myround(mod[tabinmod,2]*100,2), ")", sep = "")))}
# Condense interaction rows into one row
interact_tab <- matrix(rep(NA,2*ncol(tab)),nrow = 2)
rownames(interact_tab) <- rep("Direct interlock wtd. DV $\\times$ Years after base year",2)
interact_idx <- seq(grep(rownames(tab),pattern = "[:]")[1],grep(rownames(tab),pattern = "[:]")[length(grep(rownames(tab),pattern = "[:]"))], by = 2) 
for(i in interact_idx){
  interact_tab[1:2,which(interact_idx==i)] <- tab[c(i,i+1),which(interact_idx==i)]
}
tab <- rbind(rbind(tab[1:(ncol(tab)*2),],interact_tab),tab[c((max(grep(rownames(tab),pattern = "[:]"))+1):nrow(tab)),])
# Fix row names
rownames(tab) <- c(rbind(vars[-grep(vars,pattern = "[:]")[-1]], "")); 
rownames(tab)[seq(1,nrow(tab),by=2)] <- 
  c("CSOs, Interlock wtd.","CDP reporting, Interlock wtd.","Climate coalitions, Interlock wtd.","Climate lobbying, Interlock wtd.",
    "Direct interlock wtd. DV $\\times$ Years after base year",
    "Structural eqv. wtd. DV","Num. interlocks", "Eig. Centrality", "Num. opp. coalitions","Employees","Revenue") 
# Add N
N <- c(sum(summary(mod1a_time)$df[1:2]),sum(summary(mod2a_time)$df[1:2]),sum(summary(mod3a_time)$df[1:2]),sum(summary(mod4a_time)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(10,11); addtorow$command <- c(paste0(' '),paste0(' '))
ptable(cbind(rownames(tab), tab), "output/time_saturated.tex")

# clean up
rems <- ls()[ls() %in% c("dat","cur","pdat","pcur","myround","ptable","dat_lag","cur_lag") == FALSE]
rm(list = rems)

# 11. Calculate percentage of within-industry interlocks -------------------------------------

# calculate interlock summary stats for pg 8.
mean(sinh(dat$num_interlocks[dat$year == "2019"]), na.rm = TRUE) # avg number of interlocks
mean(sinh(dat$num_interlocks[dat$year == "2019"][order(dat$revt[dat$year == "2019"], decreasing = TRUE)][1:500]))
mean(dat$num_interlocks[dat$year == "2019"][order(dat$revt[dat$year == "2019"], decreasing = TRUE)][1:500] > 0)

# calculate interlock summary stats for pg 21.
mean(cur$cso_exists_ilwtd[cur$cso_exists_ilwtd > 0] == asinh(1))
mean(dat$cdp_report_ilwtd[dat$cdp_report_ilwtd > 0] == asinh(1))
mean(dat$numsupcoal_ilwtd[dat$numsupcoal_ilwtd > 0] == asinh(1))
mean(dat$lob_ilwtd[dat$lob_ilwtd > 0] == asinh(1))

mean(sinh(cur$cso_exists_ilwtd))
mean(sinh(dat$cdp_report_ilwtd))
mean(sinh(dat$numsupcoal_ilwtd))
mean(sinh(dat$lob_ilwtd))

mean(sinh(cur$cso_exists_ilwtd[cur$year > 2014]))
mean(sinh(dat$cdp_report_ilwtd[dat$year > 2014]))
mean(sinh(dat$numsupcoal_ilwtd[dat$year > 2014]))
mean(sinh(dat$lob_ilwtd[dat$year > 2014]))

# firms that do an action and then stop on pg 15.
cso_start_stop <- c()
for(i in 1:length(unique(dat$gvkey))){ # 
  # i <- 1729
  curvec <- dat$cso_exists[dat$gvkey == unique(dat$gvkey)[i]]
  if(sum(curvec, na.rm = TRUE) == 0) cso_start_stop[i] <- NA
  if(sum(curvec, na.rm = TRUE) > 0){
    curvec_ones <- c(1:length(curvec))[curvec == 1]
    if(0 %in% curvec[curvec_ones[1]:length(curvec)]) cso_start_stop[i] <- 1
    if(0 %in% curvec[curvec_ones[1]:length(curvec)] == FALSE) cso_start_stop[i] <- 0
  }
}
mean(cso_start_stop, na.rm = TRUE)

cdp_start_stop <- c()
for(i in 1:length(unique(dat$gvkey))){ # 
  # i <- 1729
  curvec <- dat$cdp_report[dat$gvkey == unique(dat$gvkey)[i]]
  if(sum(curvec, na.rm = TRUE) == 0) cdp_start_stop[i] <- NA
  if(sum(curvec, na.rm = TRUE) > 0){
    curvec_ones <- c(1:length(curvec))[curvec == 1]
    curvec_ones <- curvec_ones[is.na(curvec_ones)==FALSE]
    if(0 %in% curvec[curvec_ones[1]:length(curvec)]) cdp_start_stop[i] <- 1
    if(0 %in% curvec[curvec_ones[1]:length(curvec)] == FALSE) cdp_start_stop[i] <- 0
  }
}
mean(cdp_start_stop, na.rm = TRUE)

# descriptives for ``Descriptive findings'' subsection.
table(cur$cso_exists)
table(cur$cso_exists)/sum(table(cur$cso_exists))
aggregate(I(100*(cso_exists == 1)) ~ year, data = cur, FUN = mean, na.rm = TRUE)
aggregate(I(100*(cso_exists == 1)) ~ I(revt > quantile(cur$revt, .8)), data = cur, FUN = mean, na.rm = TRUE)
cur_lag$new_cso <- cur_lag$cso_exists - cur_lag$cso_exists_lag
cur_lag$new_cso[cur_lag$new_cso==-1] <- 0
table(cur_lag$new_cso)
table(cur_lag$new_cso[cur_lag$cso_exists_ilwtd > 0])[2]/table(cur_lag$new_cso)[2]
mean(cur_lag$new_cso[cur_lag$cso_exists_ilwtd > 0], na.rm = TRUE)/mean(cur_lag$new_cso[cur_lag$cso_exists_ilwtd == 0], na.rm = TRUE)

table(dat$cdp_report)
table(dat$cdp_report)/sum(table(dat$cdp_report))
aggregate(I(100*(cdp_report == 1)) ~ year, data = dat, FUN = mean, na.rm = TRUE)
aggregate(I(100*(cdp_report == 1)) ~ I(revt > quantile(revt, .8)), data = dat, FUN = mean, na.rm = TRUE)
dat_lag$new_cdp <- dat_lag$cdp_report - dat_lag$cdp_report_lag
dat_lag$new_cdp[dat_lag$new_cdp==-1] <- 0
table(dat_lag$new_cdp)
table(dat_lag$new_cdp[dat_lag$cdp_report_ilwtd > 0])[2]/table(dat_lag$new_cdp)[2]
mean(dat_lag$new_cdp[dat_lag$cdp_report_ilwtd > 0], na.rm = TRUE)/mean(dat_lag$new_cdp[dat_lag$cdp_report_ilwtd == 0], na.rm = TRUE)

table(I(dat$numsupcoal>0))
table(dat$numsupcoal > 0)/sum(table(dat$numsupcoal))
aggregate(I(100*(numsupcoal > 1)) ~ year, data = dat, FUN = mean, na.rm = TRUE)
aggregate(I(100*(numsupcoal > 1)) ~ I(revt > quantile(revt, .8)), data = dat, FUN = mean, na.rm = TRUE)
dat_lag$new_col <- I(dat_lag$numsupcoal > 0) - I(dat_lag$numsupcoal_lag > 0)
dat_lag$new_col[dat_lag$new_col==-1] <- 0
table(dat_lag$new_col)
table(dat_lag$new_col[dat_lag$numsupcoal_ilwtd > 0])[2]/table(dat_lag$new_col)[2]
mean(dat_lag$new_col[dat_lag$numsupcoal_ilwtd > 0], na.rm = TRUE)/mean(dat_lag$new_col[dat_lag$numsupcoal_ilwtd == 0], na.rm = TRUE)

table(dat$lob)
table(dat$lob > 0)/sum(table(dat$lob))
aggregate(I(100*(lob == 1)) ~ year, data = dat, FUN = mean, na.rm = TRUE)
aggregate(I(100*(lob == 1)) ~ I(revt > quantile(revt, .8)), data = dat, FUN = mean, na.rm = TRUE)
dat_lag$new_lob <- I(dat_lag$lob > 0) - I(dat_lag$lob_lag > 0)
dat_lag$new_lob[dat_lag$new_lob==-1] <- 0
table(dat_lag$new_lob)
table(dat_lag$new_lob[dat_lag$lob_ilwtd > 0])[2]/table(dat_lag$new_lob)[2]
mean(dat_lag$new_lob[dat_lag$lob_ilwtd > 0], na.rm = TRUE)/mean(dat_lag$new_lob[dat_lag$lob_ilwtd == 0], na.rm = TRUE)

## calculate within-industry interlocks
load("working_files/interlocks_df_temp")
# interlocks <- read_csv("data/Interlocks/interlocks_df_temp.csv") %>%
#   arrange(receiver_ticker, source_ticker, year) %>%
#   filter(receiver_ticker %in% dat$ticker & source_ticker %in% dat$ticker)
interlocks <- interlocks_df
names(interlocks)[2:3] <- c("source_gvkey","receiver_gvkey")

names(interlocks)
# Identify source and receiver NAICS
interlocks <- left_join(interlocks,
                        dat %>%
                          rename(source_gvkey = gvkey,
                                 source_naics = naics) %>%
                          select(source_gvkey, source_naics) %>%
                          distinct())

interlocks <- left_join(interlocks,
                        dat %>%
                          rename(receiver_gvkey = gvkey,
                                 receiver_naics = naics) %>%
                          select(receiver_gvkey, receiver_naics) %>%
                          distinct())
interlocks <- interlocks[is.na(interlocks$source_naics) == FALSE & is.na(interlocks$receiver_naics) == FALSE,]
interlocks <- interlocks[interlocks$year %in% c(1995:2019),]

interlocks <- interlocks %>%
  mutate(in_industry_naics_4 = if_else(nchar(source_naics)>=4 & nchar(receiver_naics)>=4 & 
                                         str_extract(source_naics,pattern = "^[0-9]{4}")==str_extract(receiver_naics,pattern = "^[0-9]{4}"),
                                       n_interlocks,0),
         in_industry_naics_3 = if_else(nchar(source_naics)>=3 & nchar(receiver_naics)>=3 & 
                                         str_extract(source_naics,pattern = "^[0-9]{3}")==str_extract(receiver_naics,pattern = "^[0-9]{3}"),
                                       n_interlocks,0),
         in_industry_naics_2 = if_else(nchar(source_naics)>=2 & nchar(receiver_naics)>=2 & 
                                         str_extract(source_naics,pattern = "^[0-9]{2}")==str_extract(receiver_naics,pattern = "^[0-9]{2}"),
                                       n_interlocks,0),
         in_industry_naics_6 = if_else(nchar(source_naics)>=6 & nchar(receiver_naics)>=6 & 
                                         str_extract(source_naics,pattern = "^[0-9]{6}")==str_extract(receiver_naics,pattern = "^[0-9]{6}"),
                                       n_interlocks,0))

cat(paste0("% of interlocks in same 2-digit NAICS: ",round(sum(interlocks$in_industry_naics_2)/sum(interlocks$n_interlocks)*100,1),"%"))
cat(paste0("% of interlocks in same 3-digit NAICS: ",round(sum(interlocks$in_industry_naics_3)/sum(interlocks$n_interlocks)*100,1),"%"))
cat(paste0("% of interlocks in same 4-digit NAICS: ",round(sum(interlocks$in_industry_naics_4)/sum(interlocks$n_interlocks)*100,1),"%"))
cat(paste0("% of interlocks in same 6-digit NAICS: ",round(sum(interlocks$in_industry_naics_6)/sum(interlocks$n_interlocks)*100,1),"%"))

# 12.  Main models not dropping no revt/no emp firms -------------------------------------
dat <- read_csv("analysis_data_all_firms.csv")

dat$ind_year <- paste(substr(dat$naics,1,2), dat$year, sep = "")
dat$ind <- substr(dat$naics,1,3)
dat$cso_exists_sewtd <- dat$cso_exists_sewtd*100000
dat$cdp_report_sewtd <- dat$cdp_report_sewtd*100000
dat$numsupcoal_sewtd <- dat$numsupcoal_sewtd*100000
dat$numoppcoal_sewtd <- dat$numoppcoal_sewtd*100000
dat$lob_sewtd <- dat$lob_sewtd*100000
dat$evc <- dat$evc/.001
dat$flpr_envrisk <- dat$flpr_envrisk/max(dat$flpr_envrisk, na.rm = TRUE)
dat$numsupcoal_d <- as.numeric(dat$numsupcoal > 0)
dat$numoppcoal_d <- as.numeric(dat$numoppcoal > 0)

dat$num_interlocks <- asinh(dat$num_interlocks)
dat$cso_exists_ilwtd <- asinh(dat$cso_exists_ilwtd)
dat$cdp_report_ilwtd <- asinh(dat$cdp_report_ilwtd)
dat$numsupcoal_ilwtd <- asinh(dat$numsupcoal_ilwtd)
dat$lob_ilwtd <- asinh(dat$lob_ilwtd)
dat$numoppcoal_ilwtd <- asinh(dat$numoppcoal_ilwtd)

dat$num_indinterlocks <- asinh(dat$num_indinterlocks)
dat$cso_exists_indilwtd <- asinh(dat$cso_exists_indilwtd)
dat$cdp_report_indilwtd <- asinh(dat$cdp_report_indilwtd)
dat$numsupcoal_indilwtd <- asinh(dat$numsupcoal_indilwtd)
dat$lob_indilwtd <- asinh(dat$lob_indilwtd)
dat$numoppcoal_indilwtd <- asinh(dat$numoppcoal_indilwtd)

# Create pdata.frame and demean by firm 
pdat <- pdata.frame(dat,index = c("gvkey","year"))
varlist <- colnames(pdat)[-which(colnames(pdat) %in% c("gvkey","conml","year","cik","cusip","naics","emissions_direct_intensity","ind_year","ind"))]
for(v in varlist){
  pdat[,v] <- Within(pdat[,v],effect = "individual",na.rm = T)
}

# Subset curtailed data frame
cur <- dat[dat$year %in% 2004:2021,]
# Create pdata.frame and demean by firm 
pcur <- pdata.frame(cur,index = c("gvkey","year"))
varlist <- colnames(pdat)[-which(colnames(pdat) %in% c("gvkey","conml","year","cik","cusip","naics","emissions_direct_intensity","ind_year","ind"))]
for(v in varlist){
  pcur[,v] <- Within(pcur[,v],effect = "individual",na.rm = T)
}

# models
mod1 <- lm(cso_exists ~ cso_exists_ilwtd, data = cur)
mod1aa <- lm(cso_exists ~ cso_exists_ilwtd + cso_exists_sewtd + num_interlocks + evc + numoppcoal + as.factor(ind) + as.factor(year), data = cur)
mod1a <- lm(cso_exists ~ cso_exists_ilwtd + cso_exists_sewtd + num_interlocks + evc + numoppcoal + as.factor(ind_year), data = cur)
mod1b <- lm(cso_exists ~ cso_exists_ilwtd + cso_exists_sewtd + num_interlocks + evc + numoppcoal + as.factor(ind_year), data = pcur)

mod2 <- lm(cdp_report ~ cdp_report_ilwtd, data = dat)
mod2aa <- lm(cdp_report ~ cdp_report_ilwtd + cdp_report_sewtd + num_interlocks + evc + numoppcoal + as.factor(ind) + as.factor(year), data = dat)
mod2a <- lm(cdp_report ~ cdp_report_ilwtd + cdp_report_sewtd + num_interlocks + evc + numoppcoal + as.factor(ind_year), data = dat)
mod2b <- lm(cdp_report ~ cdp_report_ilwtd + cdp_report_sewtd + num_interlocks + evc + numoppcoal + as.factor(ind_year), data = pdat)

mod3 <- lm(I(numsupcoal>0) ~ numsupcoal_ilwtd, data = dat)
mod3aa <- lm(I(numsupcoal>0) ~ numsupcoal_ilwtd + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + as.factor(ind) + as.factor(year), data = dat)
mod3a <- lm(I(numsupcoal>0) ~ numsupcoal_ilwtd + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + as.factor(ind_year), data = dat)
mod3b <- lm(numsupcoal_d ~ numsupcoal_ilwtd + numsupcoal_sewtd + num_interlocks + evc + numoppcoal + as.factor(ind_year), data = pdat)

mod4 <- lm(lob ~ lob_ilwtd, data = dat)
mod4aa <- lm(lob ~ lob_ilwtd + lob_sewtd + num_interlocks + evc + numoppcoal + as.factor(ind) + as.factor(year), data = dat)
mod4a <- lm(lob ~ lob_ilwtd + lob_sewtd + num_interlocks + evc + numoppcoal + as.factor(ind_year), data = dat)
mod4b <- lm(lob ~ lob_ilwtd + lob_sewtd + num_interlocks + evc + numoppcoal + as.factor(ind_year), data = pdat)

# GLMs for Table 2
mod1_glm <- glm(cso_exists ~ cso_exists_ilwtd, data = cur, family = "binomial")
mod1aa_glm <- glm(cso_exists ~ cso_exists_ilwtd + cso_exists_sewtd + num_interlocks + numoppcoal, data = cur, family = "binomial")
mod2_glm <- glm(cdp_report ~ cdp_report_ilwtd, data = dat, family = "binomial")
mod2aa_glm <- glm(cdp_report ~ cdp_report_ilwtd + cdp_report_sewtd + num_interlocks + evc + numoppcoal, data = dat, family = "binomial")
mod3_glm <- glm(I(numsupcoal>0) ~ numsupcoal_ilwtd, data = dat, maxit = 100, family = "binomial")
mod3aa_glm <- glm(I(numsupcoal>0) ~ numsupcoal_ilwtd + numsupcoal_sewtd + num_interlocks + evc + numoppcoal, data = dat, maxit = 100, family = "binomial")
mod4_glm <- glm(lob ~ lob_ilwtd, data = dat, family = "binomial")
mod4aa_glm <- glm(lob ~ lob_ilwtd + lob_sewtd + num_interlocks + evc + numoppcoal, data = dat, family = "binomial")
mod5_glm <- glm(I(numoppcoal>0) ~ numoppcoal_ilwtd, data = dat, family = "binomial")
mod5aa_glm <- glm(I(numoppcoal>0) ~ numoppcoal_ilwtd + numoppcoal_sewtd + num_interlocks + evc + numsupcoal, data = dat, family = "binomial")

coefs_CSO <- c(summary(mod1)$coef[2,1],summary(mod1aa)$coef[2,1],summary(mod1a)$coef[2,1],summary(mod1b)$coef[2,1],summary(mod1_glm)$coef[2,1],summary(mod1aa_glm)$coef[2,1])
se_CSO <- c(summary(mod1)$coef[2,2],summary(mod1aa)$coef[2,2],summary(mod1a)$coef[2,2],summary(mod1b)$coef[2,2],summary(mod1_glm)$coef[2,2],summary(mod1aa_glm)$coef[2,2])
coefs_CSO[1:4] <- 100*coefs_CSO[1:4]; se_CSO[1:4] <- 100*se_CSO[1:4]
se_CSO <- paste("(", myround(se_CSO,2), ")", sep = "")
p_CSO <- c(summary(mod1)$coef[2,4],summary(mod1aa)$coef[2,4],summary(mod1a)$coef[2,4],summary(mod1b)$coef[2,4],summary(mod1_glm)$coef[2,4],summary(mod1aa_glm)$coef[2,4])
stars <- rep("", length(coefs_CSO)); stars[p_CSO < .1] <- "^{+}"; stars[p_CSO < .05] <- "^{*}"; stars[p_CSO < .01] <- "^{**}"; stars[p_CSO < .001] <- "^{***}";  
coefs_CSO <- paste(myround(coefs_CSO,2), stars, sep = "")

coefs_CDP <- c(summary(mod2)$coef[2,1],summary(mod2aa)$coef[2,1],summary(mod2a)$coef[2,1],summary(mod2b)$coef[2,1],summary(mod2_glm)$coef[2,1],summary(mod2aa_glm)$coef[2,1])
se_CDP <- c(summary(mod2)$coef[2,2],summary(mod2aa)$coef[2,2],summary(mod2a)$coef[2,2],summary(mod2b)$coef[2,2],summary(mod2_glm)$coef[2,2],summary(mod2aa_glm)$coef[2,2])
coefs_CDP[1:4] <- 100*coefs_CDP[1:4]; se_CDP[1:4] <- 100*se_CDP[1:4]
se_CDP <- paste("(", myround(se_CDP,2), ")", sep = "")
p_CDP <- c(summary(mod2)$coef[2,4],summary(mod2aa)$coef[2,4],summary(mod2a)$coef[2,4],summary(mod2b)$coef[2,4],summary(mod2_glm)$coef[2,4],summary(mod2aa_glm)$coef[2,4])
stars <- rep("", length(coefs_CDP)); stars[p_CDP < .1] <- "^{+}"; stars[p_CDP < .05] <- "^{*}"; stars[p_CDP < .01] <- "^{**}"; stars[p_CDP < .001] <- "^{***}";  
coefs_CDP <- paste(myround(coefs_CDP,2), stars, sep = "")

coefs_COL <- c(summary(mod3)$coef[2,1],summary(mod3aa)$coef[2,1],summary(mod3a)$coef[2,1],summary(mod3b)$coef[2,1],summary(mod3_glm)$coef[2,1],summary(mod3aa_glm)$coef[2,1])
se_COL <- c(summary(mod3)$coef[2,2],summary(mod3aa)$coef[2,2],summary(mod3a)$coef[2,2],summary(mod3b)$coef[2,2],summary(mod3_glm)$coef[2,2],summary(mod3aa_glm)$coef[2,2])
coefs_COL[1:4] <- 100*coefs_COL[1:4]; se_COL[1:4] <- 100*se_COL[1:4]
se_COL <- paste("(", myround(se_COL,2), ")", sep = "")
p_COL <- c(summary(mod3)$coef[2,4],summary(mod3aa)$coef[2,4],summary(mod3a)$coef[2,4],summary(mod3b)$coef[2,4],summary(mod3_glm)$coef[2,4],summary(mod3aa_glm)$coef[2,4])
stars <- rep("", length(coefs_COL)); stars[p_COL < .1] <- "^{+}"; stars[p_COL < .05] <- "^{*}"; stars[p_COL < .01] <- "^{**}"; stars[p_COL < .001] <- "^{***}";  
coefs_COL <- paste(myround(coefs_COL,2), stars, sep = "")

coefs_LOB <- c(summary(mod4)$coef[2,1],summary(mod4aa)$coef[2,1],summary(mod4a)$coef[2,1],summary(mod4b)$coef[2,1],summary(mod4_glm)$coef[2,1],summary(mod4aa_glm)$coef[2,1])
se_LOB <- c(summary(mod4)$coef[2,2],summary(mod4aa)$coef[2,2],summary(mod4a)$coef[2,2],summary(mod4b)$coef[2,2],summary(mod4_glm)$coef[2,2],summary(mod4aa_glm)$coef[2,2])
coefs_LOB[1:4] <- 100*coefs_LOB[1:4]; se_LOB[1:4] <- 100*se_LOB[1:4]
se_LOB <- paste("(", myround(se_LOB,2), ")", sep = "")
p_LOB <- c(summary(mod4)$coef[2,4],summary(mod4aa)$coef[2,4],summary(mod4a)$coef[2,4],summary(mod4b)$coef[2,4],summary(mod4_glm)$coef[2,4],summary(mod4aa_glm)$coef[2,4])
stars <- rep("", length(coefs_LOB)); stars[p_LOB < .1] <- "^{+}"; stars[p_LOB < .05] <- "^{*}"; stars[p_LOB < .01] <- "^{**}"; stars[p_LOB < .001] <- "^{***}";  
coefs_LOB <- paste(myround(coefs_LOB,2), stars, sep = "")

addtorow <- list(); addtorow$pos <- list(1); addtorow$command <- c("")
ptable(rbind(c("CSO, Interlock wtd.", coefs_CSO),c("", se_CSO)), "output/tab2_CSO_keepnorevtfirms.tex")
addtorow <- list(); addtorow$pos <- list(1); addtorow$command <- c("")
ptable(rbind(c("CDP reporting, Interlock wtd.", coefs_CDP),c("", se_CDP)), "output/tab2_CDP_keepnorevtfirms.tex")
addtorow <- list(); addtorow$pos <- list(1); addtorow$command <- c("")
ptable(rbind(c("Climate coalitions, Interlock wtd.", coefs_COL),c("", se_COL)), "output/tab2_COL_keepnorevtfirms.tex")
addtorow <- list(); addtorow$pos <- list(1); addtorow$command <- c("")
ptable(rbind(c("Climate lobbying, Interlock wtd.", coefs_LOB),c("", se_LOB)), "output/tab2_LOB_keepnorevtfirms.tex")

# clean up
rems <- ls()[ls() %in% c("dat","cur","pdat","pcur","myround","ptable","dat_lag","cur_lag") == FALSE]
rm(list = rems)


